module cldwat2m_micro

!---------------------------------------------------------------------------------
! Purpose:
!   CAM Interface for microphysics
!
! Author: Andrew Gettelman, Hugh Morrison.
! Contributions from: Xiaohong Liu and Steve Ghan
! December 2005-May 2010
! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
!                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)         
! for questions contact Hugh Morrison, Andrew Gettelman
! e-mail: morrison@ucar.edu, andrew@ucar.edu
!---------------------------------------------------------------------------------
! modification for sub-columns, HM, (orig 8/11/10)
! This is done using the logical 'sub_column' set to .true. = subcolumns
!---------------------------------------------------------------------------------

  use shr_kind_mod,  only: r8=>shr_kind_r8
  use spmd_utils,    only: masterproc
  use ppgrid,        only: pcols, pver, pverp
  use physconst,     only: gravit, rair, tmelt, cpair, rh2o, rhoh2o
  use physconst,     only: latvap, latice
  use wv_saturation,  only: polysvp, epsqs 
  use cam_history,    only: addfld, add_default, phys_decomp, outfld
  use cam_history_support, only : fillvalue
  use cam_logfile,    only: iulog
  use phys_control,   only: phys_getopts
  use cldwat2m_macro, only: rhmini

  implicit none
  private
  save

  logical, public :: liu_in = .true.   ! True = Liu et al 2007 Ice nucleation 
                                       ! False = cooper fixed ice nucleation (MG2008)

   public :: ini_micro, mmicro_pcond,gamma

!constants remaped
  real(r8), private::  g              !gravity
  real(r8), private::  r              !Dry air Gas constant
  real(r8), private::  rv             !water vapor gas contstant
  real(r8), private::  cpp            !specific heat of dry air
  real(r8), private::  rhow           !density of liquid water
  real(r8), private::  xxlv           ! latent heat of vaporization
  real(r8), private::  xlf            !latent heat of freezing
  real(r8), private::  xxls           !latent heat of sublimation

real(r8), private:: rhosn  ! bulk density snow
real(r8), private:: rhoi   ! bulk density ice

real(r8), private:: ac,bc,as,bs,ai,bi,ar,br  !fall speed parameters 
real(r8), private:: ci,di    !ice mass-diameter relation parameters
real(r8), private:: cs,ds    !snow mass-diameter relation parameters
real(r8), private:: cr,dr    !drop mass-diameter relation parameters
real(r8), private:: f1s,f2s  !ventilation param for snow
real(r8), private:: Eii      !collection efficiency aggregation of ice
real(r8), private:: Ecc      !collection efficiency
real(r8), private:: Ecr      !collection efficiency cloud droplets/rain
real(r8), private:: f1r,f2r  !ventilation param for rain
real(r8), private:: DCS      !autoconversion size threshold
real(r8), private:: qsmall   !min mixing ratio 
real(r8), private:: bimm,aimm !immersion freezing
real(r8), private:: rhosu     !typical 850mn air density
real(r8), private:: mi0       ! new crystal mass
real(r8), private:: rin       ! radius of contact nuclei
real(r8), private:: qcvar     ! 1/relative variance of sub-grid qc
real(r8), private:: pi       ! pi

! Additional constants to help speed up code

real(r8), private:: cons1
real(r8), private:: cons2
real(r8), private:: cons3
real(r8), private:: cons4
real(r8), private:: cons5
real(r8), private:: cons6
real(r8), private:: cons7
real(r8), private:: cons8
real(r8), private:: cons9
real(r8), private:: cons10
real(r8), private:: cons11
real(r8), private:: cons12
real(r8), private:: cons13
real(r8), private:: cons14
real(r8), private:: cons15
real(r8), private:: cons16
real(r8), private:: cons17
real(r8), private:: cons18
real(r8), private:: cons19
real(r8), private:: cons20
real(r8), private:: cons21
real(r8), private:: cons22
real(r8), private:: cons23
real(r8), private:: cons24
real(r8), private:: cons25
real(r8), private:: cons27
real(r8), private:: cons28

real(r8), private:: lammini
real(r8), private:: lammaxi
real(r8), private:: lamminr
real(r8), private:: lammaxr
real(r8), private:: lammins
real(r8), private:: lammaxs

! parameters for snow/rain fraction for convective clouds
real(r8), private:: tmax_fsnow ! max temperature for transition to convective snow
real(r8), private:: tmin_fsnow ! min temperature for transition to convective snow

!needed for findsp
real(r8), private:: tt0       ! Freezing temperature

real(r8), private:: csmin,csmax,minrefl,mindbz

contains

!===============================================================================

subroutine ini_micro

!----------------------------------------------------------------------- 
! 
! Purpose: 
! initialize constants for the morrison microphysics
! called from stratiform.F90
! 
! Author: Andrew Gettelman Dec 2005
! 
!-----------------------------------------------------------------------

   integer k

   integer l,m, iaer
   real(r8) surften       ! surface tension of water w/respect to air (N/m)
   real(r8) arg

   character(len=16) :: eddy_scheme = ' '
   logical           :: history_microphysics     ! output variables for microphysics diagnostics package

   ! Query the PBL eddy scheme
   call phys_getopts(eddy_scheme_out              = eddy_scheme,           &
                     history_microphysics_out     = history_microphysics   )

   ! diagnostic precip
   call addfld ('QRAIN   ','kg/kg   ',pver, 'A','Diagnostic grid-mean rain mixing ratio'         ,phys_decomp)	
   call addfld ('QSNOW   ','kg/kg   ',pver, 'A','Diagnostic grid-mean snow mixing ratio'         ,phys_decomp)	
   call addfld ('NRAIN   ','m-3     ',pver, 'A','Diagnostic grid-mean rain number conc'         ,phys_decomp)	
   call addfld ('NSNOW   ','m-3     ',pver, 'A','Diagnostic grid-mean snow number conc'         ,phys_decomp)	

   ! size of precip
   call addfld ('RERCLD   ','m      ',pver, 'A','Diagnostic effective radius of Liquid Cloud and Rain' ,phys_decomp)

   call addfld ('DSNOW   ','m       ',pver, 'A','Diagnostic grid-mean snow diameter'         ,phys_decomp)	

   ! diagnostic radar reflectivity, cloud-averaged 
   call addfld ('REFL  ','DBz  ',pver, 'A','94 GHz radar reflectivity'       ,phys_decomp)
   call addfld ('AREFL  ','DBz  ',pver, 'A','Average 94 GHz radar reflectivity'       ,phys_decomp)
   call addfld ('FREFL  ','fraction  ',pver, 'A','Fractional occurance of radar reflectivity'       ,phys_decomp)
   
   call addfld ('CSRFL  ','DBz  ',pver, 'A','94 GHz radar reflectivity (CloudSat thresholds)'       ,phys_decomp)
   call addfld ('ACSRFL  ','DBz  ',pver, 'A','Average 94 GHz radar reflectivity (CloudSat thresholds)'       ,phys_decomp)
   call addfld ('FCSRFL  ','fraction  ',pver, 'A','Fractional occurance of radar reflectivity (CloudSat thresholds)' &
	,phys_decomp)
 
   call addfld ('AREFLZ ','mm^6/m^3 ',pver, 'A','Average 94 GHz radar reflectivity'       ,phys_decomp)

   ! Aerosol information
    call addfld ('NCAL    ','#/m3   ',pver, 'A','Number Concentation Activated for Liquid',phys_decomp)
    call addfld ('NCAI    ','#/m3   ',pver, 'A','Number Concentation Activated for Ice',phys_decomp)


   ! Average rain and snow mixing ratio (Q), number (N) and diameter (D), with frequency
   call addfld ('AQRAIN   ','kg/kg   ',pver, 'A','Average rain mixing ratio'         ,phys_decomp)	
   call addfld ('AQSNOW   ','kg/kg   ',pver, 'A','Average snow mixing ratio'         ,phys_decomp)	
   call addfld ('ANRAIN   ','m-3     ',pver, 'A','Average rain number conc'         ,phys_decomp)	
   call addfld ('ANSNOW   ','m-3     ',pver, 'A','Average snow number conc'         ,phys_decomp)
   call addfld ('ADRAIN   ','Micron  ',pver, 'A','Average rain effective Diameter'         ,phys_decomp)	
   call addfld ('ADSNOW   ','Micron  ',pver, 'A','Average snow effective Diameter'         ,phys_decomp)
   call addfld ('FREQR  ','fraction  ',pver, 'A','Fractional occurance of rain'       ,phys_decomp)
   call addfld ('FREQS  ','fraction  ',pver, 'A','Fractional occurance of snow'       ,phys_decomp)

   if ( history_microphysics) then
      call add_default ('AQSNOW   ', 1, ' ')
      call add_default ('FREQR    ', 1, ' ')
      call add_default ('FREQS    ', 1, ' ')
      call add_default ('AQRAIN   ', 1, ' ')
      call add_default ('AQSNOW   ', 1, ' ')
      call add_default ('ANRAIN   ', 1, ' ')
      call add_default ('ANSNOW   ', 1, ' ')
  end if

!declarations for morrison codes (transforms variable names)

   g= gravit                  !gravity
   r= rair       	      !Dry air Gas constant: note units(phys_constants are in J/K/kmol)
   rv= rh2o                   !water vapor gas contstant
   cpp = cpair                 !specific heat of dry air
   rhow = rhoh2o              !density of liquid water

! latent heats

   xxlv = latvap         ! latent heat vaporization
   xlf = latice          ! latent heat freezing
   xxls = xxlv + xlf     ! latent heat of sublimation

! parameters for snow/rain fraction for convective clouds

   tmax_fsnow = tmelt
   tmin_fsnow = tmelt-5._r8

! parameters below from Reisner et al. (1998)
! density parameters (kg/m3)

      rhosn = 250._r8    ! bulk density snow  (++ ceh)
      rhoi = 500._r8     ! bulk density ice
      rhow = 1000._r8    ! bulk density liquid


! fall speed parameters, V = aD^b
! V is in m/s

! droplets
	ac = 3.e7_r8
	bc = 2._r8

! snow
	as = 11.72_r8
	bs = 0.41_r8

! cloud ice
	ai = 700._r8
	bi = 1._r8

! rain
	ar = 841.99667_r8
	br = 0.8_r8

! particle mass-diameter relationship
! currently we assume spherical particles for cloud ice/snow
! m = cD^d

        pi= 3.1415927_r8

! cloud ice mass-diameter relationship

	ci = rhoi*pi/6._r8
	di = 3._r8

! snow mass-diameter relationship

	cs = rhosn*pi/6._r8
	ds = 3._r8

! drop mass-diameter relationship

	cr = rhow*pi/6._r8
	dr = 3._r8

! ventilation parameters for snow
! hall and prupacher

	f1s = 0.86_r8
	f2s = 0.28_r8

! collection efficiency, aggregation of cloud ice and snow

	Eii = 0.1_r8

! collection efficiency, accretion of cloud water by rain

	Ecr = 1.0_r8

! ventilation constants for rain

	f1r = 0.78_r8
	f2r = 0.32_r8

! autoconversion size threshold for cloud ice to snow (m)

       Dcs = 400.e-6_r8

! smallest mixing ratio considered in microphysics

	qsmall = 1.e-18_r8  

! immersion freezing parameters, bigg 1953

	bimm = 100._r8
	aimm = 0.66_r8

! typical air density at 850 mb

	rhosu = 85000._r8/(rair * tmelt)

! mass of new crystal due to aerosol freezing and growth (kg)

	mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)

! radius of contact nuclei aerosol (m)

        rin = 0.1e-6_r8

! 1 / relative variance of sub-grid cloud water distribution
! see morrison and gettelman, 2007, J. Climate for details

	qcvar = 2._r8

! freezing temperature
	tt0=273.15_r8

      pi=4._r8*atan(1.0_r8)

      !Range of cloudsat reflectivities (dBz) for analytic simulator
      csmin= -30._r8
      csmax= 26._r8
      mindbz = -99._r8
      !      minrefl = 10._r8**(mindbz/10._r8)
      minrefl = 1.26e-10_r8

! Define constants to help speed up code (limit calls to gamma function)

	cons1=gamma(1._r8+di)
	cons2=gamma(qcvar+2.47_r8)
	cons3=gamma(qcvar)
	cons4=gamma(1._r8+br)
	cons5=gamma(4._r8+br)
	cons6=gamma(1._r8+ds)
	cons7=gamma(1._r8+bs)     
	cons8=gamma(4._r8+bs)     
	cons9=gamma(qcvar+2._r8)     
	cons10=gamma(qcvar+1._r8)
	cons11=gamma(3._r8+bs)
	cons12=gamma(qcvar+1.15_r8)
	cons13=gamma(5._r8/2._r8+br/2._r8)
	cons14=gamma(5._r8/2._r8+bs/2._r8)
	cons15=gamma(qcvar+bc/3._r8)
	cons16=gamma(1._r8+bi)
	cons17=gamma(4._r8+bi)
	cons18=qcvar**2.47_r8
	cons19=qcvar**2
	cons20=qcvar**1.15_r8
	cons21=qcvar**(bc/3._r8)
	cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)	
	cons23=dcs**3
	cons24=dcs**2
	cons25=dcs**bs
	cons27=xxlv**2
	cons28=xxls**2

        lammaxi = 1._r8/10.e-6_r8
        lammini = 1._r8/(2._r8*dcs)
	lammaxr = 1._r8/20.e-6_r8
	lamminr = 1._r8/500.e-6_r8
	lammaxs = 1._r8/10.e-6_r8
	lammins = 1._r8/2000.e-6_r8

   return
 end subroutine ini_micro

!===============================================================================
!microphysics routine for each timestep goes here...

subroutine mmicro_pcond ( sub_column,       &
   lchnk, ncol, deltatin, tn,               &
   qn, qc, qi,                              &
   nc, ni, p, pdel, cldn,                   &
   liqcldf, icecldf,                        &
   cldo,                                    &
   rate1ord_cw2pr_st,                       &   
   naai, npccnin, rndst,nacon,              &
   tlat, qvlat,        &
   qctend, qitend, nctend, nitend, effc,    &
   effc_fn, effi, prect, preci,             &  
   nevapr, evapsnow,      &
   prain, prodsnow, cmeout, deffi, pgamrad, &
   lamcrad,qsout,dsout, &
         rflx,sflx, qrout,reff_rain,reff_snow,  &
   qcsevap,qisevap,qvres,cmeiout, &
   vtrmc,vtrmi,qcsedten,qisedten, &
   prao,prco,mnuccco,mnuccto,msacwio,psacwso,&
   bergso,bergo,melto,homoo,qcreso,prcio,praio,qireso,&
   mnuccro,pracso,meltsdt,frzrdt,mnuccdo)

!Author: Hugh Morrison, Andrew Gettelman, NCAR
! e-mail: morrison@ucar.edu, andrew@ucar.edu

   use wv_saturation, only: vqsatd_water

   ! input arguments
   logical,  intent(in) :: sub_column  ! True = configure for sub-columns  False = use w/o sub-columns (standard)
   integer,  intent(in) :: lchnk
   integer,  intent(in) :: ncol
   real(r8), intent(in) :: deltatin             ! time step (s)
   real(r8), intent(in) :: tn(pcols,pver)       ! input temperature (K)
   real(r8), intent(in) :: qn(pcols,pver)       ! input h20 vapor mixing ratio (kg/kg)

   ! note: all input cloud variables are grid-averaged
   real(r8), intent(inout) :: qc(pcols,pver)    ! cloud water mixing ratio (kg/kg)
   real(r8), intent(inout) :: qi(pcols,pver)    ! cloud ice mixing ratio (kg/kg)
   real(r8), intent(inout) :: nc(pcols,pver)    ! cloud water number conc (1/kg)
   real(r8), intent(inout) :: ni(pcols,pver)    ! cloud ice number conc (1/kg)
   real(r8), intent(in) :: p(pcols,pver)        ! air pressure (pa)
   real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across level (pa)
   real(r8), intent(in) :: cldn(pcols,pver)     ! cloud fraction
   real(r8), intent(in) :: icecldf(pcols,pver)  ! ice cloud fraction   
   real(r8), intent(in) :: liqcldf(pcols,pver)  ! liquid cloud fraction
   real(r8), intent(inout) :: cldo(pcols,pver)  ! old cloud fraction

   real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for direct cw to precip conversion 
                                                          ! used for scavenging
! Inputs for aerosol activation
   real(r8), intent(in) :: naai(pcols,pver)      ! ice nulceation number (from microp_aero_ts) 
   real(r8), intent(in) :: npccnin(pcols,pver)   ! ccn activated number tendency (from microp_aero_ts)
   real(r8), intent(in) :: rndst(pcols,pver,4)   ! radius of 4 dust bins for contact freezing (from microp_aero_ts)
   real(r8), intent(in) :: nacon(pcols,pver,4)   ! number in 4 dust bins for contact freezing  (from microp_aero_ts)

   ! output arguments

   real(r8), intent(out) :: tlat(pcols,pver)    ! latent heating rate       (W/kg)
   real(r8), intent(out) :: qvlat(pcols,pver)   ! microphysical tendency qv (1/s)
   real(r8), intent(out) :: qctend(pcols,pver)  ! microphysical tendency qc (1/s) 
   real(r8), intent(out) :: qitend(pcols,pver)  ! microphysical tendency qi (1/s)
   real(r8), intent(out) :: nctend(pcols,pver)  ! microphysical tendency nc (1/(kg*s))
   real(r8), intent(out) :: nitend(pcols,pver)  ! microphysical tendency ni (1/(kg*s))
   real(r8), intent(out) :: effc(pcols,pver)    ! droplet effective radius (micron)
   real(r8), intent(out) :: effc_fn(pcols,pver) ! droplet effective radius, assuming nc = 1.e8 kg-1
   real(r8), intent(out) :: effi(pcols,pver)    ! cloud ice effective radius (micron)
   real(r8), intent(out) :: prect(pcols)        ! surface precip rate (m/s)
   real(r8), intent(out) :: preci(pcols)        ! cloud ice/snow precip rate (m/s)
   real(r8), intent(out) :: nevapr(pcols,pver)  ! evaporation rate of rain + snow
   real(r8), intent(out) :: evapsnow(pcols,pver)! sublimation rate of snow
   real(r8), intent(out) :: prain(pcols,pver)   ! production of rain + snow
   real(r8), intent(out) :: prodsnow(pcols,pver)! production of snow
   real(r8), intent(out) :: cmeout(pcols,pver)  ! evap/sub of cloud
   real(r8), intent(out) :: deffi(pcols,pver)   ! ice effective diameter for optics (radiation)
   real(r8), intent(out) :: pgamrad(pcols,pver) ! ice gamma parameter for optics (radiation)
   real(r8), intent(out) :: lamcrad(pcols,pver) ! slope of droplet distribution for optics (radiation)
   real(r8), intent(out) :: rflx(pcols,pver+1)  ! grid-box average rain flux (kg m^-2 s^-1)
   real(r8), intent(out) :: sflx(pcols,pver+1)  ! grid-box average snow flux (kg m^-2 s^-1)
   real(r8), intent(out) :: qcsevap(pcols,pver) ! cloud water evaporation due to sedimentation
   real(r8), intent(out) :: qisevap(pcols,pver) ! cloud ice sublimation due to sublimation
   real(r8), intent(out) :: qvres(pcols,pver) ! residual condensation term to ensure RH < 100%
   real(r8), intent(out) :: cmeiout(pcols,pver) ! grid-mean cloud ice sub/dep
   real(r8), intent(out) :: vtrmc(pcols,pver) ! mass-weighted cloud water fallspeed
   real(r8), intent(out) :: vtrmi(pcols,pver) ! mass-weighted cloud ice fallspeed
   real(r8), intent(out) :: qcsedten(pcols,pver) ! qc sedimentation tendency
   real(r8), intent(out) :: qisedten(pcols,pver) ! qi sedimentation tendency
! microphysical process rates for output (mixing ratio tendencies)
   real(r8), intent(out) :: prao(pcols,pver) ! accretion of cloud by rain 
   real(r8), intent(out) :: prco(pcols,pver) ! autoconversion of cloud to rain
   real(r8), intent(out) :: mnuccco(pcols,pver) ! mixing rat tend due to immersion freezing
   real(r8), intent(out) :: mnuccto(pcols,pver) ! mixing ratio tend due to contact freezing
   real(r8), intent(out) :: msacwio(pcols,pver) ! mixing ratio tend due to H-M splintering
   real(r8), intent(out) :: psacwso(pcols,pver) ! collection of cloud water by snow
   real(r8), intent(out) :: bergso(pcols,pver) ! bergeron process on snow
   real(r8), intent(out) :: bergo(pcols,pver) ! bergeron process on cloud ice
   real(r8), intent(out) :: melto(pcols,pver) ! melting of cloud ice
   real(r8), intent(out) :: homoo(pcols,pver) ! homogeneos freezign cloud water
   real(r8), intent(out) :: qcreso(pcols,pver) ! residual cloud condensation due to removal of excess supersat
   real(r8), intent(out) :: prcio(pcols,pver) ! autoconversion of cloud ice to snow
   real(r8), intent(out) :: praio(pcols,pver) ! accretion of cloud ice by snow
   real(r8), intent(out) :: qireso(pcols,pver) ! residual ice deposition due to removal of excess supersat
   real(r8), intent(out) :: mnuccro(pcols,pver) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
   real(r8), intent(out) :: pracso (pcols,pver) ! mixing ratio tendency due to accretion of rain by snow (1/s)
   real(r8), intent(out) :: meltsdt(pcols,pver) ! latent heating rate due to melting of snow  (W/kg)
   real(r8), intent(out) :: frzrdt (pcols,pver) ! latent heating rate due to homogeneous freezing of rain (W/kg)
   real(r8), intent(out) :: mnuccdo(pcols,pver) ! mass tendency from ice nucleation

! local workspace
! all units mks unless otherwise stated

! temporary variables for sub-stepping 
	real(r8) :: t1(pcols,pver)
	real(r8) :: q1(pcols,pver)
	real(r8) :: qc1(pcols,pver)
   	real(r8) :: qi1(pcols,pver)
   	real(r8) :: nc1(pcols,pver)
   	real(r8) :: ni1(pcols,pver)
  	real(r8) :: tlat1(pcols,pver)
   	real(r8) :: qvlat1(pcols,pver)
   	real(r8) :: qctend1(pcols,pver)
   	real(r8) :: qitend1(pcols,pver)
   	real(r8) :: nctend1(pcols,pver)
   	real(r8) :: nitend1(pcols,pver)
   	real(r8) :: prect1(pcols)
   	real(r8) :: preci1(pcols)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   	real(r8) :: deltat        ! sub-time step (s)
        real(r8) :: omsm    ! number near unity for round-off issues
        real(r8) :: dto2    ! dt/2 (s)
        real(r8) :: mincld  ! minimum allowed cloud fraction
        real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg)
        real(r8) :: t(pcols,pver) ! temperature (K)
        real(r8) :: rho(pcols,pver) ! air density (kg m-3)
        real(r8) :: dv(pcols,pver)  ! diffusivity of water vapor in air
        real(r8) :: mu(pcols,pver)  ! viscocity of air
        real(r8) :: sc(pcols,pver)  ! schmidt number
        real(r8) :: kap(pcols,pver) ! thermal conductivity of air
        real(r8) :: rhof(pcols,pver) ! air density correction factor for fallspeed
        real(r8) :: cldmax(pcols,pver) ! precip fraction assuming maximum overlap
        real(r8) :: cldm(pcols,pver)   ! cloud fraction
        real(r8) :: icldm(pcols,pver)   ! ice cloud fraction
        real(r8) :: lcldm(pcols,pver)   ! liq cloud fraction
        real(r8) :: icwc(pcols)    ! in cloud water content (liquid+ice)
        real(r8) :: calpha(pcols)  ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cbeta(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cbetah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cgamma(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cgamah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: rcgama(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cmec1(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cmec2(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cmec3(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: cmec4(pcols) ! parameter for cond/evap (Zhang et al. 2003)
        real(r8) :: qtmp ! dummy qv 
        real(r8) :: dum  ! temporary dummy variable

        real(r8) :: cme(pcols,pver)  ! total (liquid+ice) cond/evap rate of cloud

        real(r8) :: cmei(pcols,pver) ! dep/sublimation rate of cloud ice
        real(r8) :: cwml(pcols,pver) ! cloud water mixing ratio
        real(r8) :: cwmi(pcols,pver) ! cloud ice mixing ratio
        real(r8) :: nnuccd(pver)   ! ice nucleation rate from deposition/cond.-freezing
        real(r8) :: mnuccd(pver)   ! mass tendency from ice nucleation
        real(r8) :: qcld              ! total cloud water
        real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
        real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
	real(r8) :: nctend_mixnuc(pcols,pver)
	real(r8) :: arg ! argument of erfc

        ! for calculation of rate1ord_cw2pr_st
        real(r8) :: qcsinksum_rate1ord(pver)   ! sum over iterations of cw to precip sink
        real(r8) :: qcsum_rate1ord(pver)    ! sum over iterations of cloud water       

        real(r8) :: alpha

        real(r8) :: dum1,dum2   !general dummy variables

        real(r8) :: npccn(pver)     ! droplet activation rate
        real(r8) :: qcic(pcols,pver) ! in-cloud cloud liquid mixing ratio
        real(r8) :: qiic(pcols,pver) ! in-cloud cloud ice mixing ratio
        real(r8) :: qniic(pcols,pver) ! in-precip snow mixing ratio
        real(r8) :: qric(pcols,pver) ! in-precip rain mixing ratio
        real(r8) :: ncic(pcols,pver) ! in-cloud droplet number conc
        real(r8) :: niic(pcols,pver) ! in-cloud cloud ice number conc
        real(r8) :: nsic(pcols,pver) ! in-precip snow number conc
        real(r8) :: nric(pcols,pver) ! in-precip rain number conc
        real(r8) :: lami(pver) ! slope of cloud ice size distr
        real(r8) :: n0i(pver) ! intercept of cloud ice size distr
        real(r8) :: lamc(pver) ! slope of cloud liquid size distr
        real(r8) :: n0c(pver) ! intercept of cloud liquid size distr
        real(r8) :: lams(pver) ! slope of snow size distr
        real(r8) :: n0s(pver) ! intercept of snow size distr
        real(r8) :: lamr(pver) ! slope of rain size distr
        real(r8) :: n0r(pver) ! intercept of rain size distr
        real(r8) :: cdist1(pver) ! size distr parameter to calculate droplet freezing
! combined size of precip & cloud drops
        real(r8) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud
        real(r8) :: arcld(pcols,pver) ! averaging control flag
        real(r8) :: Actmp  !area cross section of drops
        real(r8) :: Artmp  !area cross section of rain

        real(r8) :: pgam(pver) ! spectral width parameter of droplet size distr
        real(r8) :: lammax  ! maximum allowed slope of size distr
        real(r8) :: lammin  ! minimum allowed slope of size distr
        real(r8) :: nacnt   ! number conc of contact ice nuclei
        real(r8) :: mnuccc(pver) ! mixing ratio tendency due to freezing of cloud water
        real(r8) :: nnuccc(pver) ! number conc tendency due to freezing of cloud water

        real(r8) :: mnucct(pver) ! mixing ratio tendency due to contact freezing of cloud water
        real(r8) :: nnucct(pver) ! number conc tendency due to contact freezing of cloud water
        real(r8) :: msacwi(pver) ! mixing ratio tendency due to HM ice multiplication
        real(r8) :: nsacwi(pver) ! number conc tendency due to HM ice multiplication

        real(r8) :: prc(pver) ! qc tendency due to autoconversion of cloud droplets
        real(r8) :: nprc(pver) ! number conc tendency due to autoconversion of cloud droplets
        real(r8) :: nprc1(pver) ! qr tendency due to autoconversion of cloud droplets
        real(r8) :: nsagg(pver) ! ns tendency due to self-aggregation of snow
        real(r8) :: dc0  ! mean size droplet size distr
        real(r8) :: ds0  ! mean size snow size distr (area weighted)
        real(r8) :: eci  ! collection efficiency for riming of snow by droplets
        real(r8) :: psacws(pver) ! mixing rat tendency due to collection of droplets by snow
        real(r8) :: npsacws(pver) ! number conc tendency due to collection of droplets by snow
        real(r8) :: uni ! number-weighted cloud ice fallspeed
        real(r8) :: umi ! mass-weighted cloud ice fallspeed
        real(r8) :: uns(pver) ! number-weighted snow fallspeed
        real(r8) :: ums(pver) ! mass-weighted snow fallspeed
        real(r8) :: unr(pver) ! number-weighted rain fallspeed
        real(r8) :: umr(pver) ! mass-weighted rain fallspeed
        real(r8) :: unc ! number-weighted cloud droplet fallspeed
        real(r8) :: umc ! mass-weighted cloud droplet fallspeed
        real(r8) :: pracs(pver) ! mixing rat tendency due to collection of rain	by snow
        real(r8) :: npracs(pver) ! number conc tendency due to collection of rain by snow
        real(r8) :: mnuccr(pver) ! mixing rat tendency due to freezing of rain
        real(r8) :: nnuccr(pver) ! number conc tendency due to freezing of rain
        real(r8) :: pra(pver) ! mixing rat tendnency due to accretion of droplets by rain
        real(r8) :: npra(pver) ! nc tendnency due to accretion of droplets by rain
        real(r8) :: nragg(pver) ! nr tendency due to self-collection of rain
        real(r8) :: prci(pver) ! mixing rat tendency due to autoconversion of cloud ice to snow
        real(r8) :: nprci(pver) ! number conc tendency due to autoconversion of cloud ice to snow
        real(r8) :: prai(pver) ! mixing rat tendency due to accretion of cloud ice by snow
        real(r8) :: nprai(pver) ! number conc tendency due to accretion of cloud ice by snow
        real(r8) :: qvs ! liquid saturation vapor mixing ratio
        real(r8) :: qvi ! ice saturation vapor mixing ratio
        real(r8) :: dqsdt ! change of sat vapor mixing ratio with temperature
        real(r8) :: dqsidt ! change of ice sat vapor mixing ratio with temperature
        real(r8) :: ab ! correction factor for rain evap to account for latent heat
        real(r8) :: qclr ! water vapor mixing ratio in clear air
        real(r8) :: abi ! correction factor for snow sublimation to account for latent heat
        real(r8) :: epss ! 1/ sat relaxation timescale for snow
        real(r8) :: epsr ! 1/ sat relaxation timescale for rain
        real(r8) :: pre(pver) ! rain mixing rat tendency due to evaporation
        real(r8) :: prds(pver) ! snow mixing rat tendency due to sublimation
        real(r8) :: qce ! dummy qc for conservation check
        real(r8) :: qie ! dummy qi for conservation check
        real(r8) :: nce ! dummy nc for conservation check
        real(r8) :: nie ! dummy ni for conservation check
        real(r8) :: ratio ! parameter for conservation check
        real(r8) :: dumc(pcols,pver) ! dummy in-cloud qc
        real(r8) :: dumnc(pcols,pver) ! dummy in-cloud nc
        real(r8) :: dumi(pcols,pver) ! dummy in-cloud qi
        real(r8) :: dumni(pcols,pver) ! dummy in-cloud ni
        real(r8) :: dums(pcols,pver) ! dummy in-cloud snow mixing rat
        real(r8) :: dumns(pcols,pver) ! dummy in-cloud snow number conc
        real(r8) :: dumr(pcols,pver) ! dummy in-cloud rain mixing rat
        real(r8) :: dumnr(pcols,pver) ! dummy in-cloud rain number conc
! below are parameters for cloud water and cloud ice sedimentation calculations
        real(r8) :: fr(pver)
        real(r8) :: fnr(pver)
        real(r8) :: fc(pver)
        real(r8) :: fnc(pver)
        real(r8) :: fi(pver)
        real(r8) :: fni(pver)
        real(r8) :: fs(pver)
        real(r8) :: fns(pver)
        real(r8) :: faloutr(pver)
        real(r8) :: faloutnr(pver)
        real(r8) :: faloutc(pver)
        real(r8) :: faloutnc(pver)
        real(r8) :: falouti(pver)
        real(r8) :: faloutni(pver)
        real(r8) :: falouts(pver)
        real(r8) :: faloutns(pver)
        real(r8) :: faltndr
        real(r8) :: faltndnr
        real(r8) :: faltndc
        real(r8) :: faltndnc
        real(r8) :: faltndi
        real(r8) :: faltndni
        real(r8) :: faltnds
        real(r8) :: faltndns
        real(r8) :: faltndqie
        real(r8) :: faltndqce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        real(r8) :: relhum(pcols,pver) ! relative humidity
        real(r8) :: csigma(pcols) ! parameter for cond/evap of cloud water/ice
        real(r8) :: rgvm ! max fallspeed for all species
        real(r8) :: arn(pcols,pver) ! air density corrected rain fallspeed parameter
        real(r8) :: asn(pcols,pver) ! air density corrected snow fallspeed parameter
        real(r8) :: acn(pcols,pver) ! air density corrected cloud droplet fallspeed parameter
        real(r8) :: ain(pcols,pver) ! air density corrected cloud ice fallspeed parameter
        real(r8) :: nsubi(pver) ! evaporation of cloud ice number
        real(r8) :: nsubc(pver) ! evaporation of droplet number
        real(r8) :: nsubs(pver) ! evaporation of snow number
        real(r8) :: nsubr(pver) ! evaporation of rain number
	real(r8) :: mtime ! factor to account for droplet activation timescale
	real(r8) :: dz(pcols,pver) ! height difference across model vertical level

!fice variable
	real(r8) :: nfice(pcols,pver)

	!! add precip flux variables for sub-stepping
	real(r8) :: rflx1(pcols,pver+1)
        real(r8) :: sflx1(pcols,pver+1)

! returns from function/subroutine calls
        real(r8) :: tsp(pcols,pver)      ! saturation temp (K)
        real(r8) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
        real(r8) :: qsphy(pcols,pver)      ! saturation mixing ratio (kg/kg): hybrid rh
        real(r8) :: qs(pcols)            ! liquid-ice weighted sat mixing rat (kg/kg)
        real(r8) :: es(pcols)            ! liquid-ice weighted sat vapor press (pa)
        real(r8) :: esl(pcols,pver)      ! liquid sat vapor pressure (pa)
        real(r8) :: esi(pcols,pver)      ! ice sat vapor pressure (pa)
        real(r8) :: gammas(pcols)        ! parameter for cond/evap of cloud water

! sum of source/sink terms for diagnostic precip

        real(r8) :: qnitend(pcols,pver) ! snow mixing ratio source/sink term
        real(r8) :: nstend(pcols,pver)  ! snow number concentration source/sink term
        real(r8) :: qrtend(pcols,pver) ! rain mixing ratio source/sink term
        real(r8) :: nrtend(pcols,pver)  ! rain number concentration source/sink term
	real(r8) :: qrtot ! vertically-integrated rain mixing rat source/sink term
	real(r8) :: nrtot ! vertically-integrated rain number conc source/sink term
	real(r8) :: qstot ! vertically-integrated snow mixing rat source/sink term
	real(r8) :: nstot ! vertically-integrated snow number conc source/sink term

! new terms for Bergeron process

	real(r8) :: dumnnuc ! provisional ice nucleation rate (for calculating bergeron)
	real(r8) :: ninew  ! provisional cloud ice number conc (for calculating bergeron)
	real(r8) :: qinew ! provisional cloud ice mixing ratio (for calculating bergeron)
	real(r8) :: qvl  ! liquid sat mixing ratio   
	real(r8) :: epsi ! 1/ sat relaxation timecale for cloud ice
	real(r8) :: prd ! provisional deposition rate of cloud ice at water sat 
	real(r8) :: berg(pcols,pver) ! mixing rat tendency due to bergeron process for cloud ice
	real(r8) :: bergs(pver) ! mixing rat tendency due to bergeron process for snow

!bergeron terms
        real(r8) :: bergtsf   !bergeron timescale to remove all liquid
        real(r8) :: rhin      !modified RH for vapor deposition

! diagnostic rain/snow for output to history
! values are in-precip (local) !!!!

        real(r8), intent(out) :: qrout(pcols,pver)     ! grid-box average rain mixing ratio (kg/kg)
	real(r8), intent(out) :: reff_rain(pcols,pver) ! rain effective radius (micron)
	real(r8), intent(out) :: reff_snow(pcols,pver) ! snow effective radius (micron)
	real(r8) 	      :: drout(pcols,pver)     ! rain diameter (m)
	real(r8) :: nrout(pcols,pver) ! rain number concentration (1/m3)
	real(r8) :: nsout(pcols,pver) ! snow number concentration (1/m3)
	real(r8), intent(out) :: dsout(pcols,pver) ! snow diameter (m)
	real(r8), intent(out) :: qsout(pcols,pver) ! snow mixing ratio (kg/kg)

!averageed rain/snow for history
	real(r8) :: qrout2(pcols,pver)
	real(r8) :: qsout2(pcols,pver)
	real(r8) :: nrout2(pcols,pver)
	real(r8) :: nsout2(pcols,pver)
	real(r8) :: freqs(pcols,pver)
	real(r8) :: freqr(pcols,pver)
	real(r8) :: dumfice
	real(r8) :: drout2(pcols,pver) ! mean rain particle diameter (m)
	real(r8) :: dsout2(pcols,pver) ! mean snow particle diameter (m)

!ice nucleation, droplet activation
        real(r8) :: dum2i(pcols,pver) ! number conc of ice nuclei available (1/kg)
        real(r8) :: dum2l(pcols,pver) ! number conc of CCN (1/kg)
	real(r8) :: ncmax
	real(r8) :: nimax

!output fields for number conc
        real(r8) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3)
        real(r8) :: ncal(pcols,pver) ! output number conc of CCN (1/m3)

! loop array variables
         integer i,k,nstep,n, l
	 integer ii,kk, m

! loop variables for sub-step solution
	 integer iter,it,ltrue(pcols)

! used in contact freezing via dust particles
        real(r8)  tcnt, viscosity, mfp
        real(r8)  slip1, slip2, slip3, slip4
!        real(r8)  dfaer1, dfaer2, dfaer3, dfaer4
!        real(r8)  nacon1,nacon2,nacon3,nacon4
        real(r8)  ndfaer1, ndfaer2, ndfaer3, ndfaer4
        real(r8)  nslip1, nslip2, nslip3, nslip4

! used in ice effective radius
        real(r8)  bbi, cci, ak, iciwc, rvi

! used in Bergeron processe and water vapor deposition
        real(r8)  Tk, deles, Aprpr, Bprpr, Cice, qi0, Crate, qidep

! mean cloud fraction over the time step
        real(r8)  cldmw(pcols,pver)

! used in secondary ice production
        real(r8) ni_secp

! variabels to check for RH after rain evap

	real(r8) :: esn
	real(r8) :: qsn
	real(r8) :: ttmp


        real(r8) :: refl(pcols,pver)    ! analytic radar reflectivity        
        real(r8) :: rainrt(pcols,pver)  ! rain rate for reflectivity calculation
        real(r8) :: rainrt1(pcols,pver)
 	real(r8) :: csrfl(pcols,pver)	!cloudsat reflectivity 
        real(r8) :: arefl(pcols,pver)  !average reflectivity will zero points outside valid range
  	real(r8) :: acsrfl(pcols,pver) !cloudsat average
 	real(r8) :: frefl(pcols,pver)
  	real(r8) :: fcsrfl(pcols,pver)
 	real(r8) :: areflz(pcols,pver)  !average reflectivity in z.
	real(r8) :: tmp

        real(r8) dmc,ssmc,dstrn  ! variables for modal scheme.
  
      real(r8), parameter :: cdnl    = 0.e6_r8    ! cloud droplet number limiter

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! initialize  output fields for number conc qand ice nucleation
    ncai(1:ncol,1:pver)=0._r8 
    ncal(1:ncol,1:pver)=0._r8  

!Initialize rain size
    rercld(1:ncol,1:pver)=0._r8
    arcld(1:ncol,1:pver)=0._r8

!initialize radiation output variables
        pgamrad(1:ncol,1:pver)=0._r8 ! liquid gamma parameter for optics (radiation)
        lamcrad(1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
        deffi  (1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
!initialize radiation output variables
!initialize water vapor tendency term output
        qcsevap(1:ncol,1:pver)=0._r8 
        qisevap(1:ncol,1:pver)=0._r8 
        qvres  (1:ncol,1:pver)=0._r8 
        cmeiout (1:ncol,1:pver)=0._r8
        vtrmc (1:ncol,1:pver)=0._r8
        vtrmi (1:ncol,1:pver)=0._r8
        qcsedten (1:ncol,1:pver)=0._r8
        qisedten (1:ncol,1:pver)=0._r8    

        prao(1:ncol,1:pver)=0._r8 
        prco(1:ncol,1:pver)=0._r8 
        mnuccco(1:ncol,1:pver)=0._r8 
        mnuccto(1:ncol,1:pver)=0._r8 
        msacwio(1:ncol,1:pver)=0._r8 
        psacwso(1:ncol,1:pver)=0._r8 
        bergso(1:ncol,1:pver)=0._r8 
        bergo(1:ncol,1:pver)=0._r8 
        melto(1:ncol,1:pver)=0._r8 
        homoo(1:ncol,1:pver)=0._r8 
        qcreso(1:ncol,1:pver)=0._r8 
        prcio(1:ncol,1:pver)=0._r8 
        praio(1:ncol,1:pver)=0._r8 
        qireso(1:ncol,1:pver)=0._r8 
        mnuccro(1:ncol,1:pver)=0._r8 
        pracso (1:ncol,1:pver)=0._r8 
        meltsdt(1:ncol,1:pver)=0._r8
        frzrdt (1:ncol,1:pver)=0._r8
        mnuccdo(1:ncol,1:pver)=0._r8

! assign variable deltat for sub-stepping...
	deltat=deltatin	

! parameters for scheme

        omsm=0.99999_r8
        dto2=0.5_r8*deltat
	mincld=0.0001_r8
 
! initialize multi-level fields
        q(1:ncol,1:pver)=qn(1:ncol,1:pver)
        t(1:ncol,1:pver)=tn(1:ncol,1:pver)
	
! initialize time-varying parameters

        do k=1,pver
           do i=1,ncol
        rho(i,k)=p(i,k)/(r*t(i,k))
	dv(i,k) = 8.794E-5_r8*t(i,k)**1.81_r8/p(i,k)
	mu(i,k) = 1.496E-6_r8*t(i,k)**1.5_r8/ &
             (t(i,k)+120._r8)	
	sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
	kap(i,k) = 1.414e3_r8*1.496e-6_r8*t(i,k)** &
            1.5_r8/(t(i,k)+120._r8) 

! air density adjustment for fallspeed parameters
! includes air density correction factor to the
! power of 0.54 following Heymsfield and Bansemer 2007

	rhof(i,k)=(rhosu/rho(i,k))**0.54

        arn(i,k)=ar*rhof(i,k)
        asn(i,k)=as*rhof(i,k)
        acn(i,k)=ac*rhof(i,k)
        ain(i,k)=ai*rhof(i,k)

! get dz from dp and hydrostatic approx
! keep dz positive (define as layer k-1 - layer k)

        dz(i,k)= pdel(i,k)/(rho(i,k)*g)

        end do
        end do

! initialization
        t1(1:ncol,1:pver) = t(1:ncol,1:pver)
	q1(1:ncol,1:pver) = q(1:ncol,1:pver)
	qc1(1:ncol,1:pver) = qc(1:ncol,1:pver)
	qi1(1:ncol,1:pver) = qi(1:ncol,1:pver)
	nc1(1:ncol,1:pver) = nc(1:ncol,1:pver)
	ni1(1:ncol,1:pver) = ni(1:ncol,1:pver)

! initialize tendencies to zero
        tlat1(1:ncol,1:pver)=0._r8
	qvlat1(1:ncol,1:pver)=0._r8
	qctend1(1:ncol,1:pver)=0._r8
	qitend1(1:ncol,1:pver)=0._r8
	nctend1(1:ncol,1:pver)=0._r8
	nitend1(1:ncol,1:pver)=0._r8

! initialize precip output
	qrout(1:ncol,1:pver)=0._r8
	qsout(1:ncol,1:pver)=0._r8
	nrout(1:ncol,1:pver)=0._r8
	nsout(1:ncol,1:pver)=0._r8
        dsout(1:ncol,1:pver)=0._r8

        drout(1:ncol,1:pver)=0._r8
        !! initialize as fillvalue to avoid Floating Exceptions
	reff_rain(1:ncol,1:pver)=fillvalue
        reff_snow(1:ncol,1:pver)=fillvalue

! initialize variables for trop_mozart
	nevapr(1:ncol,1:pver) = 0._r8
	evapsnow(1:ncol,1:pver) = 0._r8
	prain(1:ncol,1:pver) = 0._r8
	prodsnow(1:ncol,1:pver) = 0._r8
	cmeout(1:ncol,1:pver) = 0._r8

! for refl calc
        rainrt1(1:ncol,1:pver) = 0._r8

! initialize precip fraction and output tendencies
        cldmax(1:ncol,1:pver)=mincld

!initialize aerosol number
!        naer2(1:ncol,1:pver,:)=0._r8
        dum2l(1:ncol,1:pver)=0._r8
        dum2i(1:ncol,1:pver)=0._r8

! initialize avg precip rate
	prect1(1:ncol)=0._r8
	preci1(1:ncol)=0._r8

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!Get humidity and saturation vapor pressures

        do k=1,pver

! find wet bulk temperature and saturation value for provisional t and q without
! condensation

        call vqsatd_water(t(1,k),p(1,k),es,qs,gammas,ncol) ! use rhw

        do i=1,ncol

	esl(i,k)=polysvp(t(i,k),0)
	esi(i,k)=polysvp(t(i,k),1)

! hm fix, make sure when above freezing that esi=esl, not active yet
	if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)

        relhum(i,k)=q(i,k)/qs(i)

! get cloud fraction, check for minimum

           cldm(i,k)=max(cldn(i,k),mincld)
           cldmw(i,k)=max(cldn(i,k),mincld)

           icldm(i,k)=max(icecldf(i,k),mincld)
           lcldm(i,k)=max(liqcldf(i,k),mincld)

! subcolumns, set cloud fraction variables to one
! if cloud water or ice is present, if not present
! set to mincld (mincld used instead of zero, to prevent
! possible division by zero errors

           if (sub_column) then

              cldm(i,k)=mincld
              cldmw(i,k)=mincld
              icldm(i,k)=mincld
              lcldm(i,k)=mincld
              
              if (qc(i,k).ge.qsmall) then
                 lcldm(i,k)=1.           
                 cldm(i,k)=1.
                 cldmw(i,k)=1.
              end if
              
              if (qi(i,k).ge.qsmall) then             
                 cldm(i,k)=1.
                 icldm(i,k)=1.
              end if
              
           end if               ! sub-columns

! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)

        nfice(i,k)=0._r8
        dumfice=qc(i,k)+qi(i,k)
        if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then
           nfice(i,k)=qi(i,k)/dumfice
        endif

	if (t(i,k).lt.tmelt - 5._r8) then
           if (liu_in) then 

! if aerosols interact with ice set number of activated ice nuclei
              dum2=naai(i,k)

           else
! cooper curve (factor of 1000 is to convert from L-1 to m-3)
              dum2=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
              dum2=min(dum,208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1
           endif

           dumnnuc=(dum2-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
           dumnnuc=max(dumnnuc,0._r8)
! get provisional ni and qi after nucleation in order to calculate
! Bergeron process below
           ninew=ni(i,k)+dumnnuc*deltat
           qinew=qi(i,k)+dumnnuc*deltat*mi0
!T>268
        else
           ninew=ni(i,k)
           qinew=qi(i,k)
        end if

! Initialize CME components

  cme(i,k) = 0._r8
  cmei(i,k)=0._r8


!-------------------------------------------------------------------
!Bergeron process

! make sure to initialize bergeron process to zero
       berg(i,k)=0._r8
       prd = 0._r8

!condensation loop.

! get in-cloud qi and ni after nucleation
       if (icldm(i,k) .gt. 0._r8) then 
          qiic(i,k)=qinew/icldm(i,k)
          niic(i,k)=ninew/icldm(i,k)
       else
          qiic(i,k)=0._r8
          niic(i,k)=0._r8
       endif

!if T < 0 C then bergeron.
           if (t(i,k).lt.273.15) then
!if ice exists
              if (qi(i,k).gt.qsmall) then

                 bergtsf = 0._r8 ! bergeron time scale (fraction of timestep)

                    qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
                    qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))

                    dqsidt =  xxls*qvi/(rv*t(i,k)**2)
                    abi = 1._r8+dqsidt*xxls/cpp

! get ice size distribution parameters

                    if (qiic(i,k).ge.qsmall) then
                       lami(k) = (cons1*ci* &
                       niic(i,k)/qiic(i,k))**(1._r8/di)
                       n0i(k) = niic(i,k)*lami(k)

! check for slope
! adjust vars
                       if (lami(k).lt.lammini) then

                          lami(k) = lammini
                          n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
                       else if (lami(k).gt.lammaxi) then
                          lami(k) = lammaxi
                          n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
                       end if
                    
                       epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k) &
                       /(lami(k)*lami(k))

!if liquid exists  
                 if (qc(i,k).gt. qsmall) then 

!begin bergeron process
!     do bergeron (vapor deposition with RHw=1)
!     code to find berg (a rate) goes here

! calculate Bergeron process

                       prd = epsi*(qvl-qvi)/abi
                    
                    else
                       prd = 0._r8
                    end if

! multiply by cloud fraction

                    prd = prd*min(icldm(i,k),lcldm(i,k))

!     transfer of existing cloud liquid to ice

                    berg(i,k)=max(0._r8,prd)

                 end if  !end liquid exists bergeron

                 if (berg(i,k).gt.0._r8) then
                    bergtsf=max(0._r8,(qc(i,k)/berg(i,k))/deltat) 

                    if(bergtsf.lt.1._r8) berg(i,k) = max(0._r8,qc(i,k)/deltat)

                 endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

                 if (bergtsf.lt.1._r8.or.icldm(i,k).gt.lcldm(i,k)) then

			if (qiic(i,k).ge.qsmall) then

! first case is for case when liquid water is present, but is completely depleted in time step, i.e., bergrsf > 0 but < 1

		if (qc(i,k).ge.qsmall) then
                       rhin  = (1.0_r8 + relhum(i,k)) / 2._r8
                       if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
                          prd = epsi*(rhin*qvl-qvi)/abi

! multiply by cloud fraction assuming liquid/ice maximum overlap
                          prd = prd*min(icldm(i,k),lcldm(i,k))

! add to cmei
                          cmei(i,k) = cmei(i,k) + (prd * (1._r8- bergtsf))

	                end if ! rhin 
	        end if ! qc > qsmall

! second case is for pure ice cloud, either no liquid, or icldm > lcldm

		if (qc(i,k).lt.qsmall.or.icldm(i,k).gt.lcldm(i,k)) then

! note: for case of no liquid, need to set liquid cloud fraction to zero
! store liquid cloud fraction in 'dum'

			if (qc(i,k).lt.qsmall) then 
			dum=0._r8 
			else
			dum=lcldm(i,k)
			end if

! set RH to grid-mean value for pure ice cloud
                       rhin = relhum(i,k)

                if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then

                       prd = epsi*(rhin*qvl-qvi)/abi

! multiply by relevant cloud fraction for pure ice cloud
! assuming maximum overlap of liquid/ice
                       prd = prd*max((icldm(i,k)-dum),0._r8)
                       cmei(i,k) = cmei(i,k) + prd

		end if ! rhin
		end if ! qc or icldm > lcldm
	end if ! qiic
	end if ! bergtsf or icldm > lcldm

!     if deposition, it should not reduce grid mean rhi below 1.0
                 if(cmei(i,k) > 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) > 1._r8 ) &
                  cmei(i,k)=min(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat)

              end if            !end ice exists loop
   !this ends temperature < 0. loop

!-------------------------------------------------------------------
           end if  ! 
!..............................................................

! evaporation should not exceed available water

           if ((-berg(i,k)).lt.-qc(i,k)/deltat) &
             berg(i,k) = max(qc(i,k)/deltat,0._r8)

!sublimation process...

            if ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall ) then

               qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
               qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
               dqsidt =  xxls*qvi/(rv*t(i,k)**2)
               abi = 1._r8+dqsidt*xxls/cpp

! get ice size distribution parameters

               lami(k) = (cons1*ci* &
               niic(i,k)/qiic(i,k))**(1._r8/di)
               n0i(k) = niic(i,k)*lami(k)

! check for slope
! adjust vars
               if (lami(k).lt.lammini) then
                  
                  lami(k) = lammini
                  n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
               else if (lami(k).gt.lammaxi) then
                  lami(k) = lammaxi
                  n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
               end if
                    
               epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k) &
               /(lami(k)*lami(k))

! modify for ice fraction below
               prd = epsi*(relhum(i,k)*qvl-qvi)/abi * icldm(i,k)
               cmei(i,k)=min(prd,0._r8)

            endif

! sublimation should not exceed available ice
           if (cmei(i,k).lt.-qi(i,k)/deltat) &
              cmei(i,k)=-qi(i,k)/deltat

! sublimation should not increase grid mean rhi above 1.0 
           if(cmei(i,k) < 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) < 1._r8 ) &
              cmei(i,k)=min(0._r8,max(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat))

! limit cmei due for roundoff error

           cmei(i,k)=cmei(i,k)*omsm

! conditional for ice nucleation 

        if (t(i,k).lt.(tmelt - 5._r8)) then 
          if ( liu_in ) then 

! using Liu et al. (2007) ice nucleation with hooks into simulated aerosol
! ice nucleation rate (dum2) has already been calculated and read in (naai)

             dum2i(i,k)=naai(i,k)
          else

! cooper curve (factor of 1000 is to convert from L-1 to m-3)
              dum2i(i,k)=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8
! put limit on number of nucleated crystals, set to number at T=-30 C
! cooper (limit to value at -35 C)
              dum2i(i,k)=min(dum2i(i,k),208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1
           endif
        else
           dum2i(i,k)=0._r8
	end if

	end do ! i loop
	end do ! k loop

 
     cldo(:ncol,:)=cldn(:ncol,:)

!! initialize sub-step precip flux variables
	do i=1,ncol
	   !! flux is zero at top interface, so these should stay as 0.
           rflx1(i,1)=0._r8
           sflx1(i,1)=0._r8
	do k=1,pver

	   ! initialize normal and sub-step precip flux variables
           rflx1(i,k+1)=0._r8
           sflx1(i,k+1)=0._r8
	end do ! i loop
	end do ! k loop
!! initialize final precip flux variables.
	do i=1,ncol
	   !! flux is zero at top interface, so these should stay as 0.
           rflx(i,1)=0._r8
           sflx(i,1)=0._r8
	do k=1,pver
	   ! initialize normal and sub-step precip flux variables
           rflx(i,k+1)=0._r8
           sflx(i,k+1)=0._r8
	end do ! i loop
	end do ! k loop	

	do i=1,ncol
	ltrue(i)=0
	do k=1,pver
! skip microphysical calculations if no cloud water

	if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall.or.cmei(i,k).ge.qsmall) ltrue(i)=1
	end do
	end do

! assign number of sub-steps to iter
! use 2 sub-steps, following tests described in MG2008
	iter = 2

! get sub-step time step
	deltat=deltat/real(iter)

! since activation/nucleation processes are fast, need to take into account
! factor mtime = mixing timescale in cloud / model time step
! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity
! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk

!        note: mtime for bulk aerosols was set to: mtime=deltat/1200._r8

        mtime=1._r8
        rate1ord_cw2pr_st(:,:)=0._r8 ! rce 2010/05/01

!!!! skip calculations if no cloud water
	do i=1,ncol
	if (ltrue(i).eq.0) then
           tlat(i,1:pver)=0._r8
           qvlat(i,1:pver)=0._r8
           qctend(i,1:pver)=0._r8
           qitend(i,1:pver)=0._r8
           qnitend(i,1:pver)=0._r8
           qrtend(i,1:pver)=0._r8
           nctend(i,1:pver)=0._r8
           nitend(i,1:pver)=0._r8
           nrtend(i,1:pver)=0._r8
           nstend(i,1:pver)=0._r8
           prect(i)=0._r8
           preci(i)=0._r8
           qniic(i,1:pver)=0._r8
           qric(i,1:pver)=0._r8
           nsic(i,1:pver)=0._r8
           nric(i,1:pver)=0._r8
           rainrt(i,1:pver)=0._r8
	goto 300
	end if

        qcsinksum_rate1ord(1:pver)=0._r8 
        qcsum_rate1ord(1:pver)=0._r8 


!!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!.....................................................................................................
	do it=1,iter

! initialize sub-step microphysical tendencies

	tlat(i,1:pver)=0._r8
	qvlat(i,1:pver)=0._r8
	qctend(i,1:pver)=0._r8
	qitend(i,1:pver)=0._r8
	qnitend(i,1:pver)=0._r8
	qrtend(i,1:pver)=0._r8
	nctend(i,1:pver)=0._r8
	nitend(i,1:pver)=0._r8
	nrtend(i,1:pver)=0._r8
	nstend(i,1:pver)=0._r8

! initialize diagnostic precipitation to zero

        qniic(i,1:pver)=0._r8
        qric(i,1:pver)=0._r8
        nsic(i,1:pver)=0._r8
        nric(i,1:pver)=0._r8
 
        rainrt(i,1:pver)=0._r8


! begin new i,k loop, calculate new cldmax after adjustment to cldm above

! initialize vertically-integrated rain and snow tendencies

	   qrtot = 0._r8
	   nrtot = 0._r8
	   qstot = 0._r8
	   nstot = 0._r8

! initialize precip at surface

	   prect(i)=0._r8
	   preci(i)=0._r8

	   do k=1,pver

! set cwml and cwmi to current qc and qi

	   cwml(i,k)=qc(i,k)
	   cwmi(i,k)=qi(i,k)

! initialize precip fallspeeds to zero

	   ums(k)=0._r8 
	   uns(k)=0._r8 
	   umr(k)=0._r8 
	   unr(k)=0._r8

! calculate precip fraction based on maximum overlap assumption

! for sub-columns cldm has already been set to 1 if cloud
! water or ice is present, so cldmax will be correctly set below
! and nothing extra needs to be done here

           if (k.eq.1) then
		cldmax(i,k)=cldm(i,k)
	   else
! if rain or snow mix ratio is smaller than
! threshold, then set cldmax to cloud fraction at current level
	   if (qric(i,k-1).ge.qsmall.or.qniic(i,k-1).ge.qsmall) then
                cldmax(i,k)=max(cldmax(i,k-1),cldm(i,k))
	   else
	   cldmax(i,k)=cldm(i,k)
	   end if
	   end if

! decrease in number concentration due to sublimation/evap
! divide by cloud fraction to get in-cloud decrease
! don't reduce Nc due to bergeron process

           if (cmei(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. cldm(i,k) > mincld) then
              nsubi(k)=cmei(i,k)/qi(i,k)*ni(i,k)/cldm(i,k)
	   else
              nsubi(k)=0._r8
           end if
           nsubc(k)=0._r8


!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%

        if (dum2i(i,k).gt.0._r8.and.t(i,k).lt.(tmelt - 5._r8).and. &
	   relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then

!if NCAI > 0. then set numice = ncai (as before)
!note: this is gridbox averaged

              nnuccd(k)=(dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
              nnuccd(k)=max(nnuccd(k),0._r8)
              nimax = dum2i(i,k)*icldm(i,k)

!Calc mass of new particles using new crystal mass...
!also this will be multiplied by mtime as nnuccd is...

              mnuccd(k) = nnuccd(k) * mi0

!  add mnuccd to cmei....
              cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime

!  limit cmei

               qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
               dqsidt =  xxls*qvi/(rv*t(i,k)**2)
               abi = 1._r8+dqsidt*xxls/cpp
              cmei(i,k)=min(cmei(i,k),(q(i,k)-qvi)/abi/deltat)

! limit for roundoff error
              cmei(i,k)=cmei(i,k)*omsm

           else
              nnuccd(k)=0._r8
	      nimax = 0._r8
              mnuccd(k) = 0._r8
           end if

!c............................................................................
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
! for microphysical process calculations
! units are kg/kg for mixing ratio, 1/kg for number conc

! limit in-cloud values to 0.005 kg/kg

           qcic(i,k)=min(cwml(i,k)/lcldm(i,k),5.e-3_r8)
           qiic(i,k)=min(cwmi(i,k)/icldm(i,k),5.e-3_r8)
           ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
           niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)

           if (qc(i,k) - berg(i,k)*deltat.lt.qsmall) then
              qcic(i,k)=0._r8
              ncic(i,k)=0._r8
              if (qc(i,k)-berg(i,k)*deltat.lt.0._r8) then
                 berg(i,k)=qc(i,k)/deltat*omsm
              end if
	   end if

	   if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.qsmall) then
	   qiic(i,k)=0._r8
	   niic(i,k)=0._r8
	   if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.0._r8) then
	   cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm
	   end if
	   end if

! add to cme output

           cmeout(i,k) = cmeout(i,k)+cmei(i,k)

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! droplet activation
! calculate potential for droplet activation if cloud water is present
! formulation from Abdul-Razzak and Ghan (2000) and Abdul-Razzak et al. (1998), AR98
! number tendency (npccnin) is read in from companion routine

! assume aerosols already activated are equal to number of existing droplets for simplicity
! multiply by cloud fraction to obtain grid-average tendency

           if (qcic(i,k).ge.qsmall) then   
              npccn(k) = max(0._r8,npccnin(i,k))  
	      dum2l(i,k)=(nc(i,k)+npccn(k)*deltat)/lcldm(i,k)
              dum2l(i,k)=max(dum2l(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3  
	      ncmax = dum2l(i,k)*lcldm(i,k)
           else
              npccn(k)=0._r8
              dum2l(i,k)=0._r8
              ncmax = 0._r8
           end if 

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! get size distribution parameters based on in-cloud cloud water/ice 
! these calculations also ensure consistency between number and mixing ratio
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!......................................................................
! cloud ice

	if (qiic(i,k).ge.qsmall) then

! add upper limit to in-cloud number concentration to prevent numerical error
	niic(i,k)=min(niic(i,k),qiic(i,k)*1.e20_r8)

	lami(k) = (cons1*ci* &
              niic(i,k)/qiic(i,k))**(1._r8/di)
	n0i(k) = niic(i,k)*lami(k)

! check for slope
! adjust vars

	if (lami(k).lt.lammini) then

	lami(k) = lammini
	n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
	niic(i,k) = n0i(k)/lami(k)
	else if (lami(k).gt.lammaxi) then
	lami(k) = lammaxi
	n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        niic(i,k) = n0i(k)/lami(k)
	end if

	else
	lami(k) = 0._r8
	n0i(k) = 0._r8
	end if

	if (qcic(i,k).ge.qsmall) then

! add upper limit to in-cloud number concentration to prevent numerical error
	ncic(i,k)=min(ncic(i,k),qcic(i,k)*1.e20_r8)

        ncic(i,k)=max(ncic(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm  

! get pgam from fit to observations of martin et al. 1994

        pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
        pgam(k)=1._r8/(pgam(k)**2)-1._r8
        pgam(k)=max(pgam(k),2._r8)
        pgam(k)=min(pgam(k),15._r8)

! calculate lamc

	lamc(k) = (pi/6._r8*rhow*ncic(i,k)*gamma(pgam(k)+4._r8)/ &
                 (qcic(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)

! lammin, 50 micron diameter max mean size

	lammin = (pgam(k)+1._r8)/50.e-6_r8
	lammax = (pgam(k)+1._r8)/2.e-6_r8

	if (lamc(k).lt.lammin) then
	lamc(k) = lammin
	ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
                gamma(pgam(k)+1._r8)/ &
               (pi*rhow*gamma(pgam(k)+4._r8))
	else if (lamc(k).gt.lammax) then
	lamc(k) = lammax
	ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
                gamma(pgam(k)+1._r8)/ &
               (pi*rhow*gamma(pgam(k)+4._r8))
	end if

! parameter to calculate droplet freezing

        cdist1(k) = ncic(i,k)/gamma(pgam(k)+1._r8) 

	else
	lamc(k) = 0._r8
        cdist1(k) = 0._r8
	end if

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! begin micropysical process calculations 
!.................................................................
! autoconversion of cloud liquid water to rain
! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
! minimum qc of 1 x 10^-8 prevents floating point error

	if (qcic(i,k).ge.1.e-8_r8) then

! nprc is increase in rain number conc due to autoconversion
! nprc1 is decrease in cloud droplet conc due to autoconversion

! assume exponential sub-grid distribution of qc, resulting in additional
! factor related to qcvar below

 ! hm switch for sub-columns, don't include sub-grid qc
           if (sub_column) then
 
              prc(k) = 1350._r8*qcic(i,k)**2.47_r8* &
              (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
              nprc(k) = prc(k)/(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
              nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
 
           else

              prc(k) = cons2/(cons3*cons18)*1350._r8*qcic(i,k)**2.47_r8* &
              (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
              nprc(k) = prc(k)/cons22
              nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
 
           end if               ! sub-column switch

        else
           prc(k)=0._r8
           nprc(k)=0._r8
	   nprc1(k)=0._r8
	end if
 
! add autoconversion to precip from above to get provisional rain mixing ratio
! and number concentration (qric and nric)

! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)

	dum=0.45_r8
	dum1=0.45_r8

	if (k.eq.1) then
	qric(i,k)=prc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
	nric(i,k)=nprc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
	else
	if (qric(i,k-1).ge.qsmall) then
	dum=umr(k-1)
	dum1=unr(k-1)
	end if

! no autoconversion of rain number if rain/snow falling from above
! this assumes that new drizzle drops formed by autoconversion are rapidly collected
! by the existing rain/snow particles from above

	if (qric(i,k-1).ge.1.e-9_r8.or.qniic(i,k-1).ge.1.e-9_r8) then
	nprc(k)=0._r8
	end if

        qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*((pra(k-1)+prc(k))*lcldm(i,k)+(pre(k-1)-pracs(k-1)-mnuccr(k-1))*cldmax(i,k))))&
                   /(dum*rho(i,k)*cldmax(i,k))
	nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*(nprc(k)*lcldm(i,k)+(nsubr(k-1)-npracs(k-1)-nnuccr(k-1)+nragg(k-1))*cldmax(i,k))))&
                   /(dum1*rho(i,k)*cldmax(i,k))

	end if

!.......................................................................
! Autoconversion of cloud ice to snow
! similar to Ferrier (1994)

	if (t(i,k).le.273.15_r8.and.qiic(i,k).ge.qsmall) then

! note: assumes autoconversion timescale of 180 sec

           nprci(k) = n0i(k)/(lami(k)*180._r8)*exp(-lami(k)*dcs)

           prci(k) = pi*rhoi*n0i(k)/(6._r8*180._r8)* &
              (cons23/lami(k)+3._r8*cons24/lami(k)**2+ &
          6._r8*dcs/lami(k)**3+6._r8/lami(k)**4)*exp(-lami(k)*dcs)          

        else
              prci(k)=0._r8
              nprci(k)=0._r8
	end if

! add autoconversion to flux from level above to get provisional snow mixing ratio
! and number concentration (qniic and nsic)

	dum=(asn(i,k)*cons25)
	dum1=(asn(i,k)*cons25)

	if (k.eq.1) then
           qniic(i,k)=prci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
           nsic(i,k)=nprci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
	else
	if (qniic(i,k-1).ge.qsmall) then
	dum=ums(k-1)
	dum1=uns(k-1)
	end if

        qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*((prci(k)+prai(k-1)+psacws(k-1)+bergs(k-1))*icldm(i,k)+(prds(k-1)+ &
         pracs(k-1)+mnuccr(k-1))*cldmax(i,k))))&
                    /(dum*rho(i,k)*cldmax(i,k))

	nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+nnuccr(k-1))*cldmax(i,k))))&
                   /(dum1*rho(i,k)*cldmax(i,k))

	end if

! if precip mix ratio is zero so should number concentration

	if (qniic(i,k).lt.qsmall) then
	qniic(i,k)=0._r8
	nsic(i,k)=0._r8
	end if

	if (qric(i,k).lt.qsmall) then
	qric(i,k)=0._r8
	nric(i,k)=0._r8
	end if

! make sure number concentration is a positive number to avoid 
! taking root of negative later

	nric(i,k)=max(nric(i,k),0._r8)
	nsic(i,k)=max(nsic(i,k),0._r8)

!.......................................................................
! get size distribution parameters for precip
!......................................................................
! rain

	if (qric(i,k).ge.qsmall) then
	lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
	n0r(k) = nric(i,k)*lamr(k)

! check for slope
! adjust vars

	if (lamr(k).lt.lamminr) then

	lamr(k) = lamminr

	n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
	nric(i,k) = n0r(k)/lamr(k)
	else if (lamr(k).gt.lammaxr) then
	lamr(k) = lammaxr
	n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
	nric(i,k) = n0r(k)/lamr(k)
	end if

! provisional rain number and mass weighted mean fallspeed (m/s)

        unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
        umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))

	else
	lamr(k) = 0._r8
	n0r(k) = 0._r8
	umr(k) = 0._r8
	unr(k) = 0._r8
	end if

!......................................................................
! snow

	if (qniic(i,k).ge.qsmall) then
	lams(k) = (cons6*cs*nsic(i,k)/ &
            qniic(i,k))**(1._r8/ds)
	n0s(k) = nsic(i,k)*lams(k)

! check for slope
! adjust vars

	if (lams(k).lt.lammins) then
	lams(k) = lammins
	n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
	nsic(i,k) = n0s(k)/lams(k)

	else if (lams(k).gt.lammaxs) then
	lams(k) = lammaxs
	n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
	nsic(i,k) = n0s(k)/lams(k)
	end if

! provisional snow number and mass weighted mean fallspeed (m/s)

        ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
        uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))

	else
	lams(k) = 0._r8
	n0s(k) = 0._r8
	ums(k) = 0._r8
	uns(k) = 0._r8
	end if

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! heterogeneous freezing of cloud water

	if (qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then

! immersion freezing (Bigg, 1953)

 
! subcolumns
           
           if (sub_column) then
 
              mnuccc(k) = &
                       pi*pi/36._r8*rhow* &
                   cdist1(k)*gamma(7._r8+pgam(k))* &
                    bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
                    lamc(k)**3/lamc(k)**3
 
              nnuccc(k) = &
                    pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
                    *bimm* &
                    (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
 
           else
 
              mnuccc(k) = cons9/(cons3*cons19)* &
                      pi*pi/36._r8*rhow* &
                  cdist1(k)*gamma(7._r8+pgam(k))* &
                   bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
                   lamc(k)**3/lamc(k)**3

              nnuccc(k) = cons10/(cons3*qcvar)* &
                  pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
                  *bimm* &
                  (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
           end if           ! sub-columns


! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
! dust size and number in 4 bins are read in from companion routine

           tcnt=(270.16_r8-t(i,k))**1.3_r8
           viscosity=1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8    ! Viscosity (kg/m/s)
           mfp=2.0_r8*viscosity/(p(i,k)  &                   ! Mean free path (m)
               *sqrt(8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i,k))))           

           nslip1=1.0_r8+(mfp/rndst(i,k,1))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,1)/mfp))))! Slip correction factor
           nslip2=1.0_r8+(mfp/rndst(i,k,2))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,2)/mfp))))
           nslip3=1.0_r8+(mfp/rndst(i,k,3))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,3)/mfp))))
           nslip4=1.0_r8+(mfp/rndst(i,k,4))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,4)/mfp))))

           ndfaer1=1.381e-23_r8*t(i,k)*nslip1/(6._r8*pi*viscosity*rndst(i,k,1))  ! aerosol diffusivity (m2/s)
           ndfaer2=1.381e-23_r8*t(i,k)*nslip2/(6._r8*pi*viscosity*rndst(i,k,2))
           ndfaer3=1.381e-23_r8*t(i,k)*nslip3/(6._r8*pi*viscosity*rndst(i,k,3))
           ndfaer4=1.381e-23_r8*t(i,k)*nslip4/(6._r8*pi*viscosity*rndst(i,k,4))


           if (sub_column) then
 
               mnucct(k) = &
                        (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
                        cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
 
               nnucct(k) = (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi*  &
                        cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
 
           else

               mnucct(k) = gamma(qcvar+4._r8/3._r8)/(cons3*qcvar**(4._r8/3._r8))*  &
                       (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
                       cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4

                nnucct(k) =  gamma(qcvar+1._r8/3._r8)/(cons3*qcvar**(1._r8/3._r8))*  &
                         (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi*  &
                       cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)

           end if      ! sub-column switch

! make sure number of droplets frozen does not exceed available ice nuclei concentration
! this prevents 'runaway' droplet freezing

        if (nnuccc(k)*lcldm(i,k).gt.nnuccd(k)) then
        dum=(nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
! scale mixing ratio of droplet freezing with limit
        mnuccc(k)=mnuccc(k)*dum
        nnuccc(k)=nnuccd(k)/lcldm(i,k)
        end if

        else
           mnuccc(k)=0._r8
           nnuccc(k)=0._r8
           mnucct(k)=0._r8
           nnucct(k)=0._r8
	end if

!.......................................................................
! snow self-aggregation from passarelli, 1978, used by reisner, 1998
! this is hard-wired for bs = 0.4 for now
! ignore self-collection of cloud ice

	if (qniic(i,k).ge.qsmall .and. t(i,k).le.273.15_r8) then
 	nsagg(k) = -1108._r8*asn(i,k)*Eii* &
             pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)*rho(i,k)** &
            ((2._r8+bs)/3._r8)*qniic(i,k)**((2._r8+bs)/3._r8)* &
            (nsic(i,k)*rho(i,k))**((4._r8-bs)/3._r8)/ &
            (4._r8*720._r8*rho(i,k))
        else
        nsagg(k)=0._r8
	end if

!.......................................................................
! accretion of cloud droplets onto snow/graupel
! here use continuous collection equation with
! simple gravitational collection kernel
! ignore collisions between droplets/cloud ice
! since minimum size ice particle for accretion is 50 - 150 micron

! ignore collision of snow with droplets above freezing

	if (qniic(i,k).ge.qsmall .and. t(i,k).le.tmelt .and. &
            qcic(i,k).ge.qsmall) then

! put in size dependent collection efficiency
! mean diameter of snow is area-weighted, since
! accretion is function of crystal geometric area
! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)

        dc0 = (pgam(k)+1._r8)/lamc(k)
        ds0 = 1._r8/lams(k)
        dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0)
        eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))

        eci = max(eci,0._r8)
        eci = min(eci,1._r8)


! no impact of sub-grid distribution of qc since psacws
! is linear in qc

	psacws(k) = pi/4._r8*asn(i,k)*qcic(i,k)*rho(i,k)* &
                  n0s(k)*Eci*cons11/ &
                  lams(k)**(bs+3._r8)	
	npsacws(k) = pi/4._r8*asn(i,k)*ncic(i,k)*rho(i,k)* &
                  n0s(k)*Eci*cons11/ &
                  lams(k)**(bs+3._r8)
        else
           psacws(k)=0._r8
           npsacws(k)=0._r8
        end if

! add secondary ice production due to accretion of droplets by snow 
! (Hallet-Mossop process) (from Cotton et al., 1986)
        if((t(i,k).lt.270.16_r8) .and. (t(i,k).ge.268.16_r8)) then
           ni_secp   = 3.5e8_r8*(270.16_r8-t(i,k))/2.0_r8*psacws(k)
           nsacwi(k) = ni_secp
           msacwi(k) = min(ni_secp*mi0,psacws(k))
        else if((t(i,k).lt.268.16_r8) .and. (t(i,k).ge.265.16_r8)) then
           ni_secp   = 3.5e8_r8*(t(i,k)-265.16_r8)/3.0_r8*psacws(k)
           nsacwi(k) = ni_secp
           msacwi(k) = min(ni_secp*mi0,psacws(k))
        else
           ni_secp   = 0.0_r8
           nsacwi(k) = 0.0_r8
           msacwi(k) = 0.0_r8
        endif
        psacws(k) = max(0.0_r8,psacws(k)-ni_secp*mi0)

!.......................................................................
! accretion of rain water by snow
! formula from ikawa and saito, 1991, used by reisner et al., 1998

	if (qric(i,k).ge.1.e-8_r8 .and. qniic(i,k).ge.1.e-8_r8 .and. & 
              t(i,k).le.273.15_r8) then

	pracs(k) = pi*pi*ecr*(((1.2_r8*umr(k)-0.95_r8*ums(k))**2+ &
                  0.08_r8*ums(k)*umr(k))**0.5_r8*rhow*rho(i,k)* &
                 n0r(k)*n0s(k)* &
                  (5._r8/(lamr(k)**6*lams(k))+ &
                  2._r8/(lamr(k)**5*lams(k)**2)+ &
                  0.5_r8/(lamr(k)**4*lams(k)**3)))

	npracs(k) = pi/2._r8*rho(i,k)*ecr*(1.7_r8*(unr(k)-uns(k))**2+ &
              0.3_r8*unr(k)*uns(k))**0.5_r8*n0r(k)*n0s(k)* &
              (1._r8/(lamr(k)**3*lams(k))+ &
              1._r8/(lamr(k)**2*lams(k)**2)+ &
              1._r8/(lamr(k)*lams(k)**3))

        else
           pracs(k)=0._r8
           npracs(k)=0._r8
	end if

!.......................................................................
! heterogeneous freezing of rain drops
! follows from Bigg (1953)

	if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then

	mnuccr(k) = 20._r8*pi*pi*rhow*nric(i,k)*bimm* &
                  (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3 &
                 /lamr(k)**3

	nnuccr(k) = pi*nric(i,k)*bimm* &
                   (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3
        else
           mnuccr(k)=0._r8
           nnuccr(k)=0._r8
	end if

!.......................................................................
! accretion of cloud liquid water by rain
! formula from Khrouditnov and Kogan (2000)
! gravitational collection kernel, droplet fall speed neglected

	if (qric(i,k).ge.qsmall .and. qcic(i,k).ge.qsmall) then

! include sub-grid distribution of cloud water

! add sub-column switch
 
           if (sub_column) then
 
              pra(k) = &
                       67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
              npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
 
           else

              pra(k) = cons12/(cons3*cons20)* &
                      67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
              npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))

           end if               ! sub-column switch

         else
           pra(k)=0._r8
           npra(k)=0._r8
	end if

!.......................................................................
! Self-collection of rain drops
! from Beheng(1994)

	if (qric(i,k).ge.qsmall) then
	nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
        else
           nragg(k)=0._r8
	end if

!.......................................................................
! Accretion of cloud ice by snow
! For this calculation, it is assumed that the Vs >> Vi
! and Ds >> Di for continuous collection

	if (qniic(i,k).ge.qsmall.and.qiic(i,k).ge.qsmall &
            .and.t(i,k).le.273.15_r8) then
	prai(k) = pi/4._r8*asn(i,k)*qiic(i,k)*rho(i,k)* &
                  n0s(k)*Eii*cons11/ &
                  lams(k)**(bs+3._r8)	
	nprai(k) = pi/4._r8*asn(i,k)*niic(i,k)* &
                  rho(i,k)*n0s(k)*Eii*cons11/ &
                  lams(k)**(bs+3._r8)
        else
        prai(k)=0._r8
        nprai(k)=0._r8
	end if

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calculate evaporation/sublimation of rain and snow
! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
! in-cloud condensation/deposition of rain and snow is neglected
! except for transfer of cloud water to snow through bergeron process

! initialize evap/sub tendncies
	pre(k)=0._r8
	prds(k)=0._r8

! evaporation of rain
! only calculate if there is some precip fraction > cloud fraction

        if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8.or.cldmax(i,k).gt.lcldm(i,k)) then

! set temporary cloud fraction to zero if cloud water + ice is very small
! this will ensure that evaporation/sublimation of precip occurs over
! entire grid cell, since min cloud fraction is specified otherwise
        if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8) then
	dum=0._r8
	else
	dum=lcldm(i,k)
	end if

! saturation vapor pressure
        esn=polysvp(t(i,k),0)
	qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)

! recalculate saturation vapor pressure for liquid and ice
	esl(i,k)=esn
	esi(i,k)=polysvp(t(i,k),1)
! hm fix, make sure when above freezing that esi=esl, not active yet
	if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)

! calculate q for out-of-cloud region
	qclr=(q(i,k)-dum*qsn)/(1._r8-dum)

	if (qric(i,k).ge.qsmall) then

        qvs=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
	dqsdt = xxlv*qvs/(rv*t(i,k)**2)
        ab = 1._r8+dqsdt*xxlv/cpp
        epsr = 2._r8*pi*n0r(k)*rho(i,k)*Dv(i,k)* &
                   (f1r/(lamr(k)*lamr(k))+ &
                    f2r*(arn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
                    sc(i,k)**(1._r8/3._r8)*cons13/ &
                (lamr(k)**(5._r8/2._r8+br/2._r8)))

	   pre(k) = epsr*(qclr-qvs)/ab

! only evaporate in out-of-cloud region
! and distribute across cldmax
           pre(k)=min(pre(k)*(cldmax(i,k)-dum),0._r8)
           pre(k)=pre(k)/cldmax(i,k)
	end if

! sublimation of snow
	if (qniic(i,k).ge.qsmall) then
        qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
        dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        abi = 1._r8+dqsidt*xxls/cpp
        epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
                   (f1s/(lams(k)*lams(k))+ &
                    f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
                    sc(i,k)**(1._r8/3._r8)*cons14/ &
               (lams(k)**(5._r8/2._r8+bs/2._r8)))
	   prds(k) = epss*(qclr-qvi)/abi	

! only sublimate in out-of-cloud region and distribute over cldmax
           prds(k)=min(prds(k)*(cldmax(i,k)-dum),0._r8)
	   prds(k)=prds(k)/cldmax(i,k)
	end if

! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
! get updated RH at end of time step based on cloud water/ice condensation/evap

	qtmp=q(i,k)-(cmei(i,k)+(pre(k)+prds(k))*cldmax(i,k))*deltat
	ttmp=t(i,k)+((pre(k)*cldmax(i,k))*xxlv+ &
                (cmei(i,k)+prds(k)*cldmax(i,k))*xxls)*deltat/cpp

!limit range of temperatures!
        ttmp=max(180._r8,min(ttmp,323._r8))

        esn=polysvp(ttmp,0)  ! use rhw to allow ice supersaturation
	qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)

! modify precip evaporation rate if q > qsat
	if (qtmp.gt.qsn) then
	if (pre(k)+prds(k).lt.-1.e-20) then
	dum1=pre(k)/(pre(k)+prds(k))
! recalculate q and t after cloud water cond but without precip evap
	qtmp=q(i,k)-(cmei(i,k))*deltat
	ttmp=t(i,k)+(cmei(i,k)*xxls)*deltat/cpp
        esn=polysvp(ttmp,0) ! use rhw to allow ice supersaturation
	qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
	dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2))
	dum=min(dum,0._r8)

! modify rates if needed, divide by cldmax to get local (in-precip) value
	pre(k)=dum*dum1/deltat/cldmax(i,k)

! do separately using RHI for prds....
        esn=polysvp(ttmp,1) ! use rhi to allow ice supersaturation
	qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)
	dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2))
	dum=min(dum,0._r8)

! modify rates if needed, divide by cldmax to get local (in-precip) value
	prds(k)=dum*(1._r8-dum1)/deltat/cldmax(i,k)
	end if
	end if

	end if

! bergeron process - evaporation of droplets and deposition onto snow

	if (qniic(i,k).ge.qsmall.and.qcic(i,k).ge.qsmall.and.t(i,k).lt.tmelt) then
        qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k))
        qvs=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(i,k))
        dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        abi = 1._r8+dqsidt*xxls/cpp
        epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
                   (f1s/(lams(k)*lams(k))+ &
                    f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
                    sc(i,k)**(1._r8/3._r8)*cons14/ &
               (lams(k)**(5._r8/2._r8+bs/2._r8)))
	bergs(k)=epss*(qvs-qvi)/abi
	else
	bergs(k)=0._r8
	end if

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! conservation to ensure no negative values of cloud water/precipitation
! in case microphysical process rates are large

! make sure and use end-of-time step values for cloud water, ice, due
! condensation/deposition

! note: for check on conservation, processes are multiplied by omsm
! to prevent problems due to round off error

! include mixing timescale  (mtime)

        qce=(qc(i,k) - berg(i,k)*deltat)
        nce=(nc(i,k)+npccn(k)*deltat*mtime)
        qie=(qi(i,k)+(cmei(i,k)+berg(i,k))*deltat)
        nie=(ni(i,k)+nnuccd(k)*deltat*mtime)

! conservation of qc

        dum = (prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+ &

              psacws(k)+bergs(k))*lcldm(i,k)*deltat

	if (dum.gt.qce) then
        ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm 

	prc(k) = prc(k)*ratio
	pra(k) = pra(k)*ratio
	mnuccc(k) = mnuccc(k)*ratio
        mnucct(k) = mnucct(k)*ratio  
        msacwi(k) = msacwi(k)*ratio  
	psacws(k) = psacws(k)*ratio	
	bergs(k) = bergs(k)*ratio	
        end if

! conservation of nc

        dum = (nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+ &
              npsacws(k)-nsubc(k))*lcldm(i,k)*deltat

	if (dum.gt.nce) then
        ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+&
            npsacws(k)-nsubc(k))*lcldm(i,k))*omsm

	nprc1(k) = nprc1(k)*ratio
	npra(k) = npra(k)*ratio
	nnuccc(k) = nnuccc(k)*ratio
        nnucct(k) = nnucct(k)*ratio  
	npsacws(k) = npsacws(k)*ratio	
	nsubc(k)=nsubc(k)*ratio
        end if

! conservation of qi

        dum = ((-mnuccc(k)-mnucct(k)-msacwi(k))*lcldm(i,k)+(prci(k)+ &
               prai(k))*icldm(i,k))*deltat

	if (dum.gt.qie) then

        ratio = (qie/deltat+(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm 
        prci(k) = prci(k)*ratio
        prai(k) = prai(k)*ratio
        end if

! conservation of ni

        dum = ((-nnucct(k)-nsacwi(k))*lcldm(i,k)+(nprci(k)+ &
               nprai(k)-nsubi(k))*icldm(i,k))*deltat

	if (dum.gt.nie) then

        ratio = (nie/deltat+(nnucct(k)+nsacwi(k))*lcldm(i,k))/ &  
                  ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
        nprci(k) = nprci(k)*ratio
        nprai(k) = nprai(k)*ratio
	nsubi(k) = nsubi(k)*ratio
        end if

! for preciptiation conservation, use logic that vertical integral 
! of tendency from current level to top of model (i.e., qrtot) cannot be negative

! conservation of rain mixing rat

	if (((prc(k)+pra(k))*lcldm(i,k)+(-mnuccr(k)+pre(k)-pracs(k))*&
               cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot.lt.0._r8) then

	if (-pre(k)+pracs(k)+mnuccr(k).ge.qsmall) then
	     
	ratio = (qrtot/(dz(i,k)*rho(i,k))+(prc(k)+pra(k))*lcldm(i,k))/&
                ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm 

        pre(k) = pre(k)*ratio
        pracs(k) = pracs(k)*ratio
        mnuccr(k) = mnuccr(k)*ratio
	end if
        end if

! conservation of nr
! for now neglect evaporation of nr
       nsubr(k)=0._r8

       if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k)&
            +nragg(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot.lt.0._r8) then

	if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k).ge.qsmall) then

	ratio = (nrtot/(dz(i,k)*rho(i,k))+nprc(k)*lcldm(i,k))/&
               ((-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k))*cldmax(i,k))*omsm

        nsubr(k) = nsubr(k)*ratio
        npracs(k) = npracs(k)*ratio
        nnuccr(k) = nnuccr(k)*ratio
        nragg(k) = nragg(k)*ratio
	end if
        end if

! conservation of snow mix ratio

	if (((bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+(pracs(k)+&
            mnuccr(k)+prds(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot.lt.0._r8) then

	if (-prds(k).ge.qsmall) then

	ratio = (qstot/(dz(i,k)*rho(i,k))+(bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+&
                (pracs(k)+mnuccr(k))*cldmax(i,k))/(-prds(k)*cldmax(i,k))*omsm

        prds(k) = prds(k)*ratio
	end if
       end if

! conservation of ns

! calculate loss of number due to sublimation
! for now neglect sublimation of ns
       nsubs(k)=0._r8

       if ((nprci(k)*icldm(i,k)+(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))*&
            dz(i,k)*rho(i,k)+nstot.lt.0._r8) then

	if (-nsubs(k)-nsagg(k).ge.qsmall) then

	ratio = (nstot/(dz(i,k)*rho(i,k))+nprci(k)*icldm(i,k)+&
                nnuccr(k)*cldmax(i,k))/((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm

        nsubs(k) = nsubs(k)*ratio
        nsagg(k) = nsagg(k)*ratio
	end if
        end if

! get tendencies due to microphysical conversion processes
! note: tendencies are multiplied by appropaiate cloud/precip 
! fraction to get grid-scale values
! note: cmei is already grid-average values

	qvlat(i,k) = qvlat(i,k)- &
       (pre(k)+prds(k))*cldmax(i,k)-cmei(i,k) 

	tlat(i,k) = tlat(i,k)+((pre(k)*cldmax(i,k)) &
           *xxlv+(prds(k)*cldmax(i,k)+cmei(i,k))*xxls+ &
           ((bergs(k)+psacws(k)+mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(mnuccr(k)+ &
               pracs(k))*cldmax(i,k)+berg(i,k))*xlf)

	qctend(i,k) = qctend(i,k)+ &
                 (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)-msacwi(k)- & 
                  psacws(k)-bergs(k))*lcldm(i,k)-berg(i,k)

	qitend(i,k) = qitend(i,k)+ &
         (mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(-prci(k)- & 
                  prai(k))*icldm(i,k)+cmei(i,k)+berg(i,k)

	qrtend(i,k) = qrtend(i,k)+ &
                 (pra(k)+prc(k))*lcldm(i,k)+(pre(k)-pracs(k)- &
                  mnuccr(k))*cldmax(i,k)

	qnitend(i,k) = qnitend(i,k)+ &
                (prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(prds(k)+ &
                   pracs(k)+mnuccr(k))*cldmax(i,k)

! add output for cmei (accumulate)
	cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)

! assign variables for trop_mozart, these are grid-average
! evaporation/sublimation is stored here as positive term

	evapsnow(i,k) = evapsnow(i,k)-prds(k)*cldmax(i,k)
	nevapr(i,k) = nevapr(i,k)-pre(k)*cldmax(i,k)

! change to make sure prain is positive: do not remove snow from
! prain used for wet deposition
        	prain(i,k) = prain(i,k)+(pra(k)+prc(k))*lcldm(i,k)+(-pracs(k)- &
                          mnuccr(k))*cldmax(i,k)
	prodsnow(i,k) = prodsnow(i,k)+(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(&
                   pracs(k)+mnuccr(k))*cldmax(i,k)

! following are used to calculate 1st order conversion rate of cloud water
!    to rain and snow (1/s), for later use in aerosol wet removal routine
! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
!    used to calculate pra, prc, ... in this routine
! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
!                      (no cloud ice or bergeron terms)
! qcsum_rate1ord     = sum over iterations{ qc used in calculation of the transfer terms }

        qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k) + (pra(k)+prc(k)+psacws(k))*lcldm(i,k) 
        qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k) 

! microphysics output, note this is grid-averaged
        prao(i,k)=prao(i,k)+pra(k)*lcldm(i,k)
        prco(i,k)=prco(i,k)+prc(k)*lcldm(i,k)
        mnuccco(i,k)=mnuccco(i,k)+mnuccc(k)*lcldm(i,k)
        mnuccto(i,k)=mnuccto(i,k)+mnucct(k)*lcldm(i,k)
        mnuccdo(i,k)=mnuccdo(i,k)+mnuccd(k)*lcldm(i,k)
        msacwio(i,k)=msacwio(i,k)+msacwi(k)*lcldm(i,k)
        psacwso(i,k)=psacwso(i,k)+psacws(k)*lcldm(i,k)
        bergso(i,k)=bergso(i,k)+bergs(k)*lcldm(i,k)
        bergo(i,k)=bergo(i,k)+berg(i,k)
        prcio(i,k)=prcio(i,k)+prci(k)*icldm(i,k)
        praio(i,k)=praio(i,k)+prai(k)*icldm(i,k)
        mnuccro(i,k)=mnuccro(i,k)+mnuccr(k)*cldmax(i,k)
        pracso (i,k)=pracso (i,k)+pracs (k)*cldmax(i,k)

! multiply activation/nucleation by mtime to account for fast timescale

	nctend(i,k) = nctend(i,k)+ npccn(k)*mtime+&
                  (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k) & 
                  -npra(k)-nprc1(k))*lcldm(i,k)      

	nitend(i,k) = nitend(i,k)+ nnuccd(k)*mtime+ & 
                  (nnucct(k)+nsacwi(k))*lcldm(i,k)+(nsubi(k)-nprci(k)- &   
                  nprai(k))*icldm(i,k)

	nstend(i,k) = nstend(i,k)+(nsubs(k)+ &
                  nsagg(k)+nnuccr(k))*cldmax(i,k)+nprci(k)*icldm(i,k)

	nrtend(i,k) = nrtend(i,k)+ &
                  nprc(k)*lcldm(i,k)+(nsubr(k)-npracs(k)-nnuccr(k) &
                   +nragg(k))*cldmax(i,k)

! make sure that nc and ni at advanced time step do not exceed
! maximum (existing N + source terms*dt), which is possible due to
! fast nucleation timescale

	if (nctend(i,k).gt.0._r8.and.nc(i,k)+nctend(i,k)*deltat.gt.ncmax) then
	nctend(i,k)=max(0._r8,(ncmax-nc(i,k))/deltat)
	end if
	if (nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax) then
	nitend(i,k)=max(0._r8,(nimax-ni(i,k))/deltat)
	end if

! get final values for precipitation q and N, based on
! flux of precip from above, source/sink term, and terminal fallspeed
! see eq. 15-16 in MG2008

! rain

        if (qric(i,k).ge.qsmall) then
	if (k.eq.1) then
	qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(k)
	nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(k)
	else
        qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(k)*rho(i,k)*cldmax(i,k))
	nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(k)*rho(i,k)*cldmax(i,k))

	end if
	else
	qric(i,k)=0._r8
	nric(i,k)=0._r8
	end if

! snow

        if (qniic(i,k).ge.qsmall) then
	if (k.eq.1) then
	qniic(i,k)=qnitend(i,k)*dz(i,k)/cldmax(i,k)/ums(k)	
	nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(k)
	else
        qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*qnitend(i,k)))/(ums(k)*rho(i,k)*cldmax(i,k))
	nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
         (rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(k)*rho(i,k)*cldmax(i,k))
	end if
	else
	qniic(i,k)=0._r8
	nsic(i,k)=0._r8
	end if

! calculate precipitation flux at surface
! divide by density of water to get units of m/s

	prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+&
                   qnitend(i,k)*dz(i,k)*rho(i,k))/rhow
	preci(i) = preci(i)+qnitend(i,k)*dz(i,k)*rho(i,k)/rhow

        ! convert rain rate from m/s to mm/hr
 
         rainrt(i,k)=qric(i,k)*rho(i,k)*umr(k)/rhow*3600._r8*1000._r8

! vertically-integrated precip source/sink terms (note: grid-averaged)

	qrtot = max(qrtot+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
	qstot = max(qstot+qnitend(i,k)*dz(i,k)*rho(i,k),0._r8)
	nrtot = max(nrtot+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
	nstot = max(nstot+nstend(i,k)*dz(i,k)*rho(i,k),0._r8)

! calculate melting and freezing of precip

! melt snow at +2 C

        if (t(i,k)+tlat(i,k)/cpp*deltat > 275.15_r8) then
           if (qstot > 0._r8) then

! make sure melting snow doesn't reduce temperature below threshold
	      dum = -xlf/cpp*qstot/(dz(i,k)*rho(i,k))
	      if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.275.15_r8) then
	       dum = (t(i,k)+tlat(i,k)/cpp*deltat-275.15_r8)*cpp/xlf
	       dum = dum/(xlf/cpp*qstot/(dz(i,k)*rho(i,k)))
	       dum = max(0._r8,dum)
	       dum = min(1._r8,dum)
	      else
               dum = 1._r8
              end if

	      qric(i,k)=qric(i,k)+dum*qniic(i,k)
	      nric(i,k)=nric(i,k)+dum*nsic(i,k)
	      qniic(i,k)=(1._r8-dum)*qniic(i,k)
	      nsic(i,k)=(1._r8-dum)*nsic(i,k)
! heating tendency 
              tmp=-xlf*dum*qstot/(dz(i,k)*rho(i,k))
              meltsdt(i,k)=meltsdt(i,k) + tmp

	      tlat(i,k)=tlat(i,k)+tmp
	      qrtot=qrtot+dum*qstot
	      nrtot=nrtot+dum*nstot
	      qstot=(1._r8-dum)*qstot
	      nstot=(1._r8-dum)*nstot
	      preci(i)=(1._r8-dum)*preci(i)
           end if
        end if

! freeze all rain at -5C for Arctic

        if (t(i,k)+tlat(i,k)/cpp*deltat < (tmelt - 5._r8)) then

           if (qrtot > 0._r8) then

! make sure freezing rain doesn't increase temperature above threshold
          dum = xlf/cpp*qrtot/(dz(i,k)*rho(i,k))
          if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.(tmelt - 5._r8)) then
           dum = -(t(i,k)+tlat(i,k)/cpp*deltat-(tmelt-5._r8))*cpp/xlf
           dum = dum/(xlf/cpp*qrtot/(dz(i,k)*rho(i,k)))
           dum = max(0._r8,dum)
           dum = min(1._r8,dum)
          else
             dum = 1._r8
          end if
            
	      qniic(i,k)=qniic(i,k)+dum*qric(i,k)
	      nsic(i,k)=nsic(i,k)+dum*nric(i,k)
	      qric(i,k)=(1._r8-dum)*qric(i,k)
	      nric(i,k)=(1._r8-dum)*nric(i,k)
! heating tendency 
	      tmp = xlf*dum*qrtot/(dz(i,k)*rho(i,k))
	      frzrdt(i,k)=frzrdt(i,k) + tmp

	      tlat(i,k)=tlat(i,k)+tmp
	      qstot=qstot+dum*qrtot
	      qrtot=(1._r8-dum)*qrtot
	      nstot=nstot+dum*nrtot
	      nrtot=(1._r8-dum)*nrtot
	      preci(i)=preci(i)+dum*(prect(i)-preci(i))
           end if
        end if

! if rain/snow mix ratio is zero so should number concentration

	if (qniic(i,k).lt.qsmall) then
	qniic(i,k)=0._r8
	nsic(i,k)=0._r8
	end if

	if (qric(i,k).lt.qsmall) then
	qric(i,k)=0._r8
	nric(i,k)=0._r8
	end if

! make sure number concentration is a positive number to avoid 
! taking root of negative

	nric(i,k)=max(nric(i,k),0._r8)
	nsic(i,k)=max(nsic(i,k),0._r8)

!.......................................................................
! get size distribution parameters for fallspeed calculations
!......................................................................
! rain

	if (qric(i,k).ge.qsmall) then
	lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
	n0r(k) = nric(i,k)*lamr(k)

! check for slope
! change lammax and lammin for rain and snow
! adjust vars

	if (lamr(k).lt.lamminr) then

	lamr(k) = lamminr

	n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
	nric(i,k) = n0r(k)/lamr(k)
	else if (lamr(k).gt.lammaxr) then
	lamr(k) = lammaxr
	n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
	nric(i,k) = n0r(k)/lamr(k)
	end if


! 'final' values of number and mass weighted mean fallspeed for rain (m/s)

        unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
        umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))

	else
	lamr(k) = 0._r8
	n0r(k) = 0._r8
	umr(k)=0._r8
	unr(k)=0._r8
	end if

!calculate mean size of combined rain and snow

        if (lamr(k).gt.0._r8) then
           Artmp = n0r(k) * pi / (2 * lamr(k)**3._r8)
        else 
           Artmp = 0._r8
        endif

        if (lamc(k).gt.0._r8) then
           Actmp = cdist1(k) * pi * gamma(pgam(k)+3._r8)/(4._r8 * lamc(k)**2._r8)
        else 
           Actmp = 0._r8
        endif

        if (Actmp.gt.0_r8.or.Artmp.gt.0) then
           rercld(i,k)=rercld(i,k) + 3._r8 *(qric(i,k) + qcic(i,k)) / (4._r8 * rhow * (Actmp + Artmp))
           arcld(i,k)=arcld(i,k)+1._r8
        endif

!......................................................................
! snow

	if (qniic(i,k).ge.qsmall) then
	lams(k) = (cons6*cs*nsic(i,k)/ &
            qniic(i,k))**(1._r8/ds)
	n0s(k) = nsic(i,k)*lams(k)

! check for slope
! adjust vars

	if (lams(k).lt.lammins) then
	lams(k) = lammins
	n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
	nsic(i,k) = n0s(k)/lams(k)

	else if (lams(k).gt.lammaxs) then
	lams(k) = lammaxs
	n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
	nsic(i,k) = n0s(k)/lams(k)
	end if

! 'final' values of number and mass weighted mean fallspeed for snow (m/s)

        ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
        uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))

	else
	lams(k) = 0._r8
	n0s(k) = 0._r8
	ums(k) = 0._r8
	uns(k) = 0._r8
	end if

!c........................................................................
! sum over sub-step for average process rates

! convert rain/snow q and N for output to history, note, 
! output is for gridbox average

	qrout(i,k)=qrout(i,k)+qric(i,k)*cldmax(i,k)
	qsout(i,k)=qsout(i,k)+qniic(i,k)*cldmax(i,k)
	nrout(i,k)=nrout(i,k)+nric(i,k)*rho(i,k)*cldmax(i,k)
	nsout(i,k)=nsout(i,k)+nsic(i,k)*rho(i,k)*cldmax(i,k)

	tlat1(i,k)=tlat1(i,k)+tlat(i,k)
	qvlat1(i,k)=qvlat1(i,k)+qvlat(i,k)
	qctend1(i,k)=qctend1(i,k)+qctend(i,k)
	qitend1(i,k)=qitend1(i,k)+qitend(i,k)
	nctend1(i,k)=nctend1(i,k)+nctend(i,k)
	nitend1(i,k)=nitend1(i,k)+nitend(i,k)

	t(i,k)=t(i,k)+tlat(i,k)*deltat/cpp
	q(i,k)=q(i,k)+qvlat(i,k)*deltat
	qc(i,k)=qc(i,k)+qctend(i,k)*deltat
	qi(i,k)=qi(i,k)+qitend(i,k)*deltat
	nc(i,k)=nc(i,k)+nctend(i,k)*deltat
	ni(i,k)=ni(i,k)+nitend(i,k)*deltat

        rainrt1(i,k)=rainrt1(i,k)+rainrt(i,k)

!divide rain radius over substeps for average
        if (arcld(i,k) .gt. 0._r8) then
           rercld(i,k)=rercld(i,k)/arcld(i,k)
        end if

!calculate precip fluxes and adding them to summing sub-stepping variables
      !! flux is zero at top interface
      	rflx(i,1)=0.0_r8
      	sflx(i,1)=0.0_r8
	
      !! calculating the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s)
      	rflx(i,k+1)=qrout(i,k)*rho(i,k)*umr(k)
	sflx(i,k+1)=qsout(i,k)*rho(i,k)*ums(k)
	
      !! add to summing sub-stepping variable
	rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1)
	sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1)	

!c........................................................................

	end do ! k loop

	prect1(i)=prect1(i)+prect(i)
	preci1(i)=preci1(i)+preci(i)

	end do ! it loop, sub-step

        do k = 1, pver               
        rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)/max(qcsum_rate1ord(k),1.0e-30_r8) 
        end do                       

 300    continue  ! continue if no cloud water
	end do ! i loop

! convert dt from sub-step back to full time step
	deltat=deltat*real(iter)

!c.............................................................................

        do i=1,ncol

! skip all calculations if no cloud water
	if (ltrue(i).eq.0) then

	do k=1,pver
! assign default values for effective radius
	effc(i,k)=10._r8
	effi(i,k)=25._r8
	effc_fn(i,k)=10._r8
        lamcrad(i,k)=0._r8
        pgamrad(i,k)=0._r8
        deffi(i,k)=0._r8
	end do
	goto 500
	end if

! initialize nstep for sedimentation sub-steps
	nstep = 1

! divide precip rate by number of sub-steps to get average over time step

	prect(i)=prect1(i)/real(iter)
	preci(i)=preci1(i)/real(iter)

       do k=1,pver

! assign variables back to start-of-timestep values before updating after sub-steps 

	t(i,k)=t1(i,k)
	q(i,k)=q1(i,k)
	qc(i,k)=qc1(i,k)
	qi(i,k)=qi1(i,k)
	nc(i,k)=nc1(i,k)
	ni(i,k)=ni1(i,k)

! divide microphysical tendencies by number of sub-steps to get average over time step

	tlat(i,k)=tlat1(i,k)/real(iter)
	qvlat(i,k)=qvlat1(i,k)/real(iter)
	qctend(i,k)=qctend1(i,k)/real(iter)
	qitend(i,k)=qitend1(i,k)/real(iter)
	nctend(i,k)=nctend1(i,k)/real(iter)
	nitend(i,k)=nitend1(i,k)/real(iter)
 
        rainrt(i,k)=rainrt1(i,k)/real(iter)

! divide by number of sub-steps to find final values
	rflx(i,k+1)=rflx1(i,k+1)/real(iter)
	sflx(i,k+1)=sflx1(i,k+1)/real(iter)

! divide output precip q and N by number of sub-steps to get average over time step

	qrout(i,k)=qrout(i,k)/real(iter)
	qsout(i,k)=qsout(i,k)/real(iter)
	nrout(i,k)=nrout(i,k)/real(iter)
	nsout(i,k)=nsout(i,k)/real(iter)

! divide trop_mozart variables by number of sub-steps to get average over time step 

	nevapr(i,k) = nevapr(i,k)/real(iter)
	evapsnow(i,k) = evapsnow(i,k)/real(iter)
	prain(i,k) = prain(i,k)/real(iter)
	prodsnow(i,k) = prodsnow(i,k)/real(iter)
	cmeout(i,k) = cmeout(i,k)/real(iter)

	cmeiout(i,k) = cmeiout(i,k)/real(iter)
        meltsdt(i,k) = meltsdt(i,k)/real(iter)
        frzrdt (i,k) = frzrdt (i,k)/real(iter)


! microphysics output
        prao(i,k)=prao(i,k)/real(iter)
        prco(i,k)=prco(i,k)/real(iter)
        mnuccco(i,k)=mnuccco(i,k)/real(iter)
        mnuccto(i,k)=mnuccto(i,k)/real(iter)
        msacwio(i,k)=msacwio(i,k)/real(iter)
        psacwso(i,k)=psacwso(i,k)/real(iter)
        bergso(i,k)=bergso(i,k)/real(iter)
        bergo(i,k)=bergo(i,k)/real(iter)
        prcio(i,k)=prcio(i,k)/real(iter)
        praio(i,k)=praio(i,k)/real(iter)

        mnuccro(i,k)=mnuccro(i,k)/real(iter)
        pracso (i,k)=pracso (i,k)/real(iter)

        mnuccdo(i,k)=mnuccdo(i,k)/real(iter)

! modify to include snow. in prain & evap (diagnostic here: for wet dep)
        nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
        prain(i,k) = prain(i,k) + prodsnow(i,k)

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! calculate sedimentation for cloud water and ice
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! update in-cloud cloud mixing ratio and number concentration 
! with microphysical tendencies to calculate sedimentation, assign to dummy vars
! note: these are in-cloud values***, hence we divide by cloud fraction

        dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
        dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
        dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
        dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)

! obtain new slope parameter to avoid possible singularity

	if (dumi(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
	dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)

	lami(k) = (cons1*ci* &
              dumni(i,k)/dumi(i,k))**(1._r8/di)
        lami(k)=max(lami(k),lammini)
        lami(k)=min(lami(k),lammaxi)
        else
           lami(k)=0._r8
        end if

	if (dumc(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
	dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
! add lower limit to in-cloud number concentration
         dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3 
         pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
         pgam(k)=1._r8/(pgam(k)**2)-1._r8
         pgam(k)=max(pgam(k),2._r8)
         pgam(k)=min(pgam(k),15._r8)

	lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
                 (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
	lammin = (pgam(k)+1._r8)/50.e-6_r8
	lammax = (pgam(k)+1._r8)/2.e-6_r8
        lamc(k)=max(lamc(k),lammin)
        lamc(k)=min(lamc(k),lammax)
        else
           lamc(k)=0._r8
        end if

! calculate number and mass weighted fall velocity for droplets
! include effects of sub-grid distribution of cloud water


	if (dumc(i,k).ge.qsmall) then
	unc = &
              acn(i,k)*gamma(1._r8+bc+pgam(k))/ &
               (lamc(k)**bc*gamma(pgam(k)+1._r8))
	umc = &
              acn(i,k)*gamma(4._r8+bc+pgam(k))/ &
             (lamc(k)**bc*gamma(pgam(k)+4._r8))
! fallspeed for output
	vtrmc(i,k)=umc
	else
	umc = 0._r8
	unc = 0._r8
	end if

! calculate number and mass weighted fall velocity for cloud ice

	if (dumi(i,k).ge.qsmall) then
	uni =  ain(i,k)*cons16/lami(k)**bi
	umi = ain(i,k)*cons17/(6._r8*lami(k)**bi)
        uni=min(uni,1.2_r8*rhof(i,k))
        umi=min(umi,1.2_r8*rhof(i,k))

! fallspeed
	vtrmi(i,k)=umi
	else
	umi = 0._r8
	uni = 0._r8
	end if

	fi(k) = g*rho(i,k)*umi
	fni(k) = g*rho(i,k)*uni
	fc(k) = g*rho(i,k)*umc
	fnc(k) = g*rho(i,k)*unc

! calculate number of split time steps to ensure courant stability criteria
! for sedimentation calculations

	rgvm = max(fi(k),fc(k),fni(k),fnc(k))
	nstep = max(int(rgvm*deltat/pdel(i,k)+1._r8),nstep)

! redefine dummy variables - sedimentation is calculated over grid-scale
! quantities to ensure conservation

        dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
        dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
        dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)

	if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
	if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8

	end do       !!! vertical loop

	do n = 1,nstep  !! loop over sub-time step to ensure stability	
	
	do k = 1,pver
	falouti(k) = fi(k)*dumi(i,k)
	faloutni(k) = fni(k)*dumni(i,k)
	faloutc(k) = fc(k)*dumc(i,k)
	faloutnc(k) = fnc(k)*dumnc(i,k)
	end do

! top of model

	k = 1
	faltndi = falouti(k)/pdel(i,k)
	faltndni = faloutni(k)/pdel(i,k)
	faltndc = faloutc(k)/pdel(i,k)
        faltndnc = faloutnc(k)/pdel(i,k)

! add fallout terms to microphysical tendencies	

	qitend(i,k) = qitend(i,k)-faltndi/nstep
	nitend(i,k) = nitend(i,k)-faltndni/nstep
	qctend(i,k) = qctend(i,k)-faltndc/nstep
        nctend(i,k) = nctend(i,k)-faltndnc/nstep

! sedimentation tendencies for output
	qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
	qisedten(i,k)=qisedten(i,k)-faltndi/nstep

	dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
	dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
	dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
        dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep

	do k = 2,pver

! for cloud liquid and ice, if cloud fraction increases with height
! then add flux from above to both vapor and cloud water of current level
! this means that flux entering clear portion of cell from above evaporates
! instantly

        dum=lcldm(i,k)/lcldm(i,k-1)
	dum=min(dum,1._r8)
        dum1=icldm(i,k)/icldm(i,k-1)
	dum1=min(dum1,1._r8)

	faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k)
	faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k)
	faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k)
	faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k)
	faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k)
	faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k)

! add fallout terms to eulerian tendencies

	qitend(i,k) = qitend(i,k)-faltndi/nstep
	nitend(i,k) = nitend(i,k)-faltndni/nstep
	qctend(i,k) = qctend(i,k)-faltndc/nstep
        nctend(i,k) = nctend(i,k)-faltndnc/nstep

! sedimentation tendencies for output
	qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
	qisedten(i,k)=qisedten(i,k)-faltndi/nstep
	
! add terms to to evap/sub of cloud water

        qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep
! for output
        qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep
        qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep
! for output
        qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep

	tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep
	tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep

	dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
	dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
	dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
        dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep

        Fni(K)=MAX(Fni(K)/pdel(i,K),Fni(K-1)/pdel(i,K-1))*pdel(i,K)
        FI(K)=MAX(FI(K)/pdel(i,K),FI(K-1)/pdel(i,K-1))*pdel(i,K)
        fnc(k)=max(fnc(k)/pdel(i,k),fnc(k-1)/pdel(i,k-1))*pdel(i,k)
        Fc(K)=MAX(Fc(K)/pdel(i,K),Fc(K-1)/pdel(i,K-1))*pdel(i,K)

	end do   !! k loop

! units below are m/s
! cloud water/ice sedimentation flux at surface 
! is added to precip flux at surface to get total precip (cloud + precip water)
! rate

	  prect(i) = prect(i)+(faloutc(pver)+falouti(pver)) &
                     /g/nstep/1000._r8	
	  preci(i) = preci(i)+(falouti(pver)) &
                     /g/nstep/1000._r8

        end do   !! nstep loop

! end sedimentation
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

! get new update for variables that includes sedimentation tendency
! note : here dum variables are grid-average, NOT in-cloud

  do k=1,pver	

        dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
        dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)

	if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
	if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
	
! calculate instantaneous processes (melting, homogeneous freezing)

        if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then
           if (dumi(i,k) > 0._r8) then

! limit so that melting does not push temperature below freezing
	      dum = -dumi(i,k)*xlf/cpp
	      if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then
	       dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf
	       dum = dum/dumi(i,k)*xlf/cpp 
	       dum = max(0._r8,dum)
	       dum = min(1._r8,dum)
	      else
	       dum = 1._r8
	      end if

              qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat

! for output
              melto(i,k)=dum*dumi(i,k)/deltat

! assume melting ice produces droplet
! mean volume radius of 8 micron

	      nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ &
               (4._r8*pi*5.12e-16_r8*rhow)

              qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat
              nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat
              tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat
           end if
        end if

! homogeneously freeze droplets at -40 C

        if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then
           if (dumc(i,k) > 0._r8) then

! limit so that freezing does not push temperature above threshold
	      dum = dumc(i,k)*xlf/cpp
	      if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then
	       dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf
	       dum = dum/dumc(i,k)*xlf/cpp
	       dum = max(0._r8,dum)
	       dum = min(1._r8,dum)
	      else
	       dum = 1._r8
	      end if

              qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat
! for output
              homoo(i,k)=dum*dumc(i,k)/deltat

! assume 25 micron mean volume radius of homogeneously frozen droplets
! consistent with size of detrained ice in stratiform.F90
	      nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* &
                 500._r8)/deltat
              qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat
              nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat
              tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat
           end if
        end if

! remove any excess over-saturation, which is possible due to non-linearity when adding 
! together all microphysical processes
! follow code similar to old CAM scheme

	qtmp=q(i,k)+qvlat(i,k)*deltat
	ttmp=t(i,k)+tlat(i,k)/cpp*deltat

        esn = polysvp(ttmp,0)  ! use rhw to allow ice supersaturation
	qsn = min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8)

	if (qtmp > qsn .and. qsn > 0) then
! expression below is approximate since there may be ice deposition
	dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat
! add to output cme
	cmeout(i,k) = cmeout(i,k)+dum
! now add to tendencies, partition between liquid and ice based on temperature
           if (ttmp > 268.15_r8) then
              dum1=0.0_r8
! now add to tendencies, partition between liquid and ice based on te
           else if (ttmp < 238.15_r8) then
              dum1=1.0_r8
           else
              dum1=(268.15_r8-ttmp)/30._r8
	   end if  

	dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
                     *qsn/(cpp*rv*ttmp**2))/deltat
	qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
! for output
        qcreso(i,k)=dum*(1._r8-dum1)
	qitend(i,k)=qitend(i,k)+dum*dum1
	qireso(i,k)=dum*dum1
	qvlat(i,k)=qvlat(i,k)-dum
! for output
        qvres(i,k)=-dum
	tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
	end if

!...............................................................................
! calculate effective radius for pass to radiation code
! if no cloud water, default value is 10 micron for droplets,
! 25 micron for cloud ice

! update cloud variables after instantaneous processes to get effective radius
! variables are in-cloud to calculate size dist parameters

        dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
        dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)

! limit in-cloud mixing ratio to reasonable value of 5 g kg-1

	dumc(i,k)=min(dumc(i,k),5.e-3_r8)
	dumi(i,k)=min(dumi(i,k),5.e-3_r8)

!...................
! cloud ice effective radius

	if (dumi(i,k).ge.qsmall) then
! add upper limit to in-cloud number concentration to prevent numerical error
	dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
	lami(k) = (cons1*ci* &
              dumni(i,k)/dumi(i,k))**(1._r8/di)

	if (lami(k).lt.lammini) then
	lami(k) = lammini
	n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
	niic(i,k) = n0i(k)/lami(k)
! adjust number conc if needed to keep mean size in reasonable range
	nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
	else if (lami(k).gt.lammaxi) then
	lami(k) = lammaxi
	n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
	niic(i,k) = n0i(k)/lami(k)
! adjust number conc if needed to keep mean size in reasonable range
	nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
	end if
	   effi(i,k) = 1.5_r8/lami(k)*1.e6_r8

        else
	   effi(i,k) = 25._r8
        end if

!...................
! cloud droplet effective radius

	if (dumc(i,k).ge.qsmall) then

! add upper limit to in-cloud number concentration to prevent numerical error
           dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! set tendency to ensure minimum droplet concentration
! after update by microphysics, except when lambda exceeds bounds on mean drop
! size or if there is no cloud water
           if (dumnc(i,k).lt.cdnl/rho(i,k)) then   
              nctend(i,k)=(cdnl/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat   
           end if
           dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
           pgam(k)=1._r8/(pgam(k)**2)-1._r8
           pgam(k)=max(pgam(k),2._r8)
           pgam(k)=min(pgam(k),15._r8)
           
           lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
                     (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
           lammin = (pgam(k)+1._r8)/50.e-6_r8
           lammax = (pgam(k)+1._r8)/2.e-6_r8
           if (lamc(k).lt.lammin) then
              lamc(k) = lammin
              ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
                   gamma(pgam(k)+1._r8)/ &
                   (pi*rhow*gamma(pgam(k)+4._r8))
! adjust number conc if needed to keep mean size in reasonable range
              nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat

           else if (lamc(k).gt.lammax) then
              lamc(k) = lammax
              ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
                   gamma(pgam(k)+1._r8)/ &
                   (pi*rhow*gamma(pgam(k)+4._r8))
! adjust number conc if needed to keep mean size in reasonable range
              nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
	end if

        effc(i,k) = &
             gamma(pgam(k)+4._r8)/ &
             gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
!assign output fields for shape here
        lamcrad(i,k)=lamc(k)
        pgamrad(i,k)=pgam(k)

        else
           effc(i,k) = 10._r8
           lamcrad(i,k)=0._r8
           pgamrad(i,k)=0._r8
        end if

! ice effective diameter for david mitchell's optics
        deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8


!!! recalculate effective radius for constant number, in order to separate
! first and second indirect effects
! assume constant number of 10^8 kg-1

	dumnc(i,k)=1.e8

	if (dumc(i,k).ge.qsmall) then
         pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
         pgam(k)=1._r8/(pgam(k)**2)-1._r8
         pgam(k)=max(pgam(k),2._r8)
         pgam(k)=min(pgam(k),15._r8)

	lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
                 (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
	lammin = (pgam(k)+1._r8)/50.e-6_r8
	lammax = (pgam(k)+1._r8)/2.e-6_r8
	if (lamc(k).lt.lammin) then
	lamc(k) = lammin
	else if (lamc(k).gt.lammax) then
	lamc(k) = lammax
	end if
	effc_fn(i,k) = &
             gamma(pgam(k)+4._r8)/ &
             gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8

        else
	effc_fn(i,k) = 10._r8
        end if


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!

	end do ! vertical k loop

 500    continue

	do k=1,pver
! if updated q (after microphysics) is zero, then ensure updated n is also zero

        if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat
        if (qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat
	end do

	end do ! i loop


! hm add rain/snow mixing ratio and number concentration as diagnostic

      call outfld('QRAIN',qrout,   pcols, lchnk)
      call outfld('QSNOW',qsout,   pcols, lchnk)
      call outfld('NRAIN',nrout,   pcols, lchnk)
      call outfld('NSNOW',nsout,   pcols, lchnk)

! add snow ouptut
      do i = 1,ncol
         do k=1,pver
            if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
               dsout(i,k)=3._r8*rhosn/917._r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
            endif
         end do
      end do
      call outfld('DSNOW',dsout,   pcols, lchnk)
 
!calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
      do i = 1,ncol
         do k=1,pver
	   !! RAIN
            if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
		reff_rain(i,k)=1.5_r8*(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)*1.e6_r8 
            endif
	    !! SNOW
            if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
		reff_snow(i,k)=1.5_r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)*1.e6_r8 

            endif
         end do
      end do

      ! analytic radar reflectivity
      ! formulas from Matthew Shupe, NOAA/CERES
      ! *****note: radar reflectivity is local (in-precip average)
      ! units of mm^6/m^3
 
      do i = 1,ncol
         do k=1,pver
            if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall) then
               dum=((qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
                  /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/cldmax(i,k)
            else
               dum=0._r8
            end if
            if (qi(i,k)+qitend(i,k)*deltat.ge.qsmall) then
               dum1=((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/cldmax(i,k)
            else 
               dum1=0._r8
            end if
         
            if (qsout(i,k).ge.qsmall) then
               dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
            end if
            
            refl(i,k)=dum+dum1
 
            ! add rain rate, but for 37 GHz formulation instead of 94 GHz
            ! formula approximated from data of Matrasov (2007)
            ! rainrt is the rain rate in mm/hr
            ! reflectivity (dum) is in DBz
 
            if (rainrt(i,k).ge.0.001) then
               dum=log10(rainrt(i,k)**6._r8)+16._r8
 
               ! convert from DBz to mm^6/m^3
 
               dum = 10._r8**(dum/10._r8)
            else
               ! don't include rain rate in R calculation for values less than 0.001 mm/hr
               dum=0.
            end if
 
            ! add to refl
 
            refl(i,k)=refl(i,k)+dum
 
            !output reflectivity in Z.
            areflz(i,k)=refl(i,k)
 
            ! convert back to DBz 
 
            if (refl(i,k).gt.minrefl) then 
               refl(i,k)=10._r8*log10(refl(i,k))
            else	
               refl(i,k)=-9999._r8
            end if
  
            !set averaging flag
            if (refl(i,k).gt.mindbz) then 
               arefl(i,k)=refl(i,k)
               frefl(i,k)=1.0_r8  
            else
               arefl(i,k)=0._r8
               areflz(i,k)=0._r8
               frefl(i,k)=0._r8
            end if
 
            ! bound cloudsat reflectivity
 
            csrfl(i,k)=min(csmax,refl(i,k))
 
            !set averaging flag
            if (csrfl(i,k).gt.csmin) then 
               acsrfl(i,k)=refl(i,k)
               fcsrfl(i,k)=1.0_r8  
            else
               acsrfl(i,k)=0._r8
               fcsrfl(i,k)=0._r8
            end if
 
         end do
      end do
       
      call outfld('REFL',refl,   pcols, lchnk)
      call outfld('AREFL',arefl,   pcols, lchnk)
      call outfld('AREFLZ',areflz,   pcols, lchnk)
      call outfld('FREFL',frefl,   pcols, lchnk)
      call outfld('CSRFL',csrfl,   pcols, lchnk)
      call outfld('ACSRFL',acsrfl,   pcols, lchnk)
      call outfld('FCSRFL',fcsrfl,   pcols, lchnk)

      call outfld('RERCLD',rercld,   pcols, lchnk)

! averaging for snow and rain number and diameter

  qrout2(:,:)=0._r8
  qsout2(:,:)=0._r8
  nrout2(:,:)=0._r8
  nsout2(:,:)=0._r8	
  drout2(:,:)=0._r8
  dsout2(:,:)=0._r8	
  freqs(:,:)=0._r8
  freqr(:,:)=0._r8
  do i = 1,ncol
     do k=1,pver
	if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
 	   qrout2(i,k)=qrout(i,k)
	   nrout2(i,k)=nrout(i,k)
	   drout2(i,k)=(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)
	   freqr(i,k)=1._r8
	endif
	if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
 	   qsout2(i,k)=qsout(i,k)
	   nsout2(i,k)=nsout(i,k)
	   dsout2(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
	   freqs(i,k)=1._r8
	endif
     end do
  end do

! output activated liquid and ice (convert from #/kg -> #/m3)
  do i = 1,ncol
     do k=1,pver
        ncai(i,k)=dum2i(i,k)*rho(i,k)
        ncal(i,k)=dum2l(i,k)*rho(i,k)
     end do
  end do
  call outfld('NCAL',ncal,    pcols,lchnk)
  call outfld('NCAI',ncai,    pcols,lchnk)

!add averaged output fields.
  call outfld('AQRAIN',qrout2,    pcols,lchnk)
  call outfld('AQSNOW',qsout2,    pcols,lchnk)
  call outfld('ANRAIN',nrout2,    pcols,lchnk)
  call outfld('ANSNOW',nsout2,    pcols,lchnk)
  call outfld('ADRAIN',drout2,    pcols,lchnk)
  call outfld('ADSNOW',dsout2,    pcols,lchnk)
  call outfld('FREQR',freqr,    pcols,lchnk)
  call outfld('FREQS',freqs,    pcols,lchnk)

!redefine fice here....
	nfice(:,:)=0._r8
	do k=1,pver
	   do i=1,ncol
        	dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        	dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
		dumfice=qsout(i,k) + qrout(i,k) + dumc(i,k) + dumi(i,k)  

                if (dumfice.gt.qsmall.and.(qsout(i,k)+dumi(i,k).gt.qsmall)) then
                  nfice(i,k)=(qsout(i,k) + dumi(i,k))/dumfice
                endif

                if (nfice(i,k).gt.1._r8) then
                   nfice(i,k)=1._r8
                endif

	   enddo
	enddo

      call outfld('FICE',nfice,   pcols, lchnk)

return
end subroutine mmicro_pcond

!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      FUNCTION GAMMA(X)

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!D    DOUBLE PRECISION FUNCTION DGAMMA(X)
!----------------------------------------------------------------------
!
! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
!   COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
!   THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
!   FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS.  COEFFICIENTS
!   FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
!   THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
!   THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
!   COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
!   MACHINE-DEPENDENT CONSTANTS.
!
!
!*******************************************************************
!*******************************************************************
!
! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
!
! BETA   - RADIX FOR THE FLOATING-POINT REPRESENTATION
! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
! XBIG   - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
!          IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
!                  GAMMA(XBIG) = BETA**MAXEXP
! XINF   - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
!          APPROXIMATELY BETA**MAXEXP
! EPS    - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
!          1.0+EPS .GT. 1.0
! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
!          1/XMININ IS MACHINE REPRESENTABLE
!
!     APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
!
!                            BETA       MAXEXP        XBIG
!
! CRAY-1         (S.P.)        2         8191        966.961
! CYBER 180/855
!   UNDER NOS    (S.P.)        2         1070        177.803
! IEEE (IBM/XT,
!   SUN, ETC.)   (S.P.)        2          128        35.040
! IEEE (IBM/XT,
!   SUN, ETC.)   (D.P.)        2         1024        171.624
! IBM 3033       (D.P.)       16           63        57.574
! VAX D-FORMAT   (D.P.)        2          127        34.844
! VAX G-FORMAT   (D.P.)        2         1023        171.489
!
!                            XINF         EPS        XMININ
!
! CRAY-1         (S.P.)   5.45E+2465   7.11E-15    1.84E-2466
! CYBER 180/855
!   UNDER NOS    (S.P.)   1.26E+322    3.55E-15    3.14E-294
! IEEE (IBM/XT,
!   SUN, ETC.)   (S.P.)   3.40E+38     1.19E-7     1.18E-38
! IEEE (IBM/XT,
!   SUN, ETC.)   (D.P.)   1.79D+308    2.22D-16    2.23D-308
! IBM 3033       (D.P.)   7.23D+75     2.22D-16    1.39D-76
! VAX D-FORMAT   (D.P.)   1.70D+38     1.39D-17    5.88D-39
! VAX G-FORMAT   (D.P.)   8.98D+307    1.11D-16    1.12D-308
!
!*******************************************************************
!*******************************************************************
!
! ERROR RETURNS
!
!  THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
!     WHEN OVERFLOW WOULD OCCUR.  THE COMPUTATION IS BELIEVED
!     TO BE FREE OF UNDERFLOW AND OVERFLOW.
!
!
!  INTRINSIC FUNCTIONS REQUIRED ARE:
!
!     INT, DBLE, EXP, LOG, REAL, SIN
!
!
! REFERENCES:  AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
!              FUNCTIONS   W. J. CODY, LECTURE NOTES IN MATHEMATICS,
!              506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
!              (ED.), SPRINGER VERLAG, BERLIN, 1976.
!
!              COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
!              SONS, NEW YORK, 1968.
!
!  LATEST MODIFICATION: OCTOBER 12, 1989
!
!  AUTHORS: W. J. CODY AND L. STOLTZ
!           APPLIED MATHEMATICS DIVISION
!           ARGONNE NATIONAL LABORATORY
!           ARGONNE, IL 60439
!
!----------------------------------------------------------------------
      INTEGER I,N
      LOGICAL PARITY

      real(r8) gamma
      REAL(r8) &
!D    DOUBLE PRECISION
         C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
         TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
      DIMENSION C(7),P(8),Q(8)
!----------------------------------------------------------------------
!  MATHEMATICAL CONSTANTS
!----------------------------------------------------------------------
      DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0_r8,0.5E0_r8,12.0E0_r8,2.0E0_r8,0.0E0_r8/, &
          SQRTPI/0.9189385332046727417803297E0_r8/, &
          PI/3.1415926535897932384626434E0_r8/
!D    DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
!D   1     SQRTPI/0.9189385332046727417803297D0/,
!D   2     PI/3.1415926535897932384626434D0/
!----------------------------------------------------------------------
!  MACHINE DEPENDENT PARAMETERS
!----------------------------------------------------------------------
      DATA XBIG,XMININ,EPS/35.040E0_r8,1.18E-38_r8,1.19E-7_r8/, &
          XINF/3.4E38_r8/
!D    DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
!D   1     XINF/1.79D308/
!----------------------------------------------------------------------
!  NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
!     APPROXIMATION OVER (1,2).
!----------------------------------------------------------------------
      DATA P/-1.71618513886549492533811E+0_r8,2.47656508055759199108314E+1_r8,&
            -3.79804256470945635097577E+2_r8,6.29331155312818442661052E+2_r8,&
            8.66966202790413211295064E+2_r8,-3.14512729688483675254357E+4_r8,&
            -3.61444134186911729807069E+4_r8,6.64561438202405440627855E+4_r8/
      DATA Q/-3.08402300119738975254353E+1_r8,3.15350626979604161529144E+2_r8,&
           -1.01515636749021914166146E+3_r8,-3.10777167157231109440444E+3_r8,&
             2.25381184209801510330112E+4_r8,4.75584627752788110767815E+3_r8,&
           -1.34659959864969306392456E+5_r8,-1.15132259675553483497211E+5_r8/
!D    DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
!D   1       -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
!D   2       8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
!D   3       -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
!D    DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
!D   1      -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
!D   2        2.25381184209801510330112D+4,4.75584627752788110767815D+3,
!D   3      -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
!----------------------------------------------------------------------
!  COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
!----------------------------------------------------------------------
      DATA C/-1.910444077728E-03_r8,8.4171387781295E-04_r8, &
          -5.952379913043012E-04_r8,7.93650793500350248E-04_r8,&
          -2.777777777777681622553E-03_r8,8.333333333333333331554247E-02_r8,&
           5.7083835261E-03_r8/
!D    DATA C/-1.910444077728D-03,8.4171387781295D-04,
!D   1     -5.952379913043012D-04,7.93650793500350248D-04,
!D   2     -2.777777777777681622553D-03,8.333333333333333331554247D-02,
!D   3      5.7083835261D-03/
!----------------------------------------------------------------------
!  STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
!----------------------------------------------------------------------
      CONV(I) = REAL(I,r8)
!D    CONV(I) = DBLE(I)
      PARITY=.FALSE.
      FACT=ONE
      N=0
      Y=X
      IF(Y.LE.ZERO)THEN
!----------------------------------------------------------------------
!  ARGUMENT IS NEGATIVE
!----------------------------------------------------------------------
        Y=-X
        Y1=AINT(Y)
        RES=Y-Y1
        IF(RES.NE.ZERO)THEN
          IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
          FACT=-PI/SIN(PI*RES)
          Y=Y+ONE
        ELSE
          RES=XINF
          GOTO 900
        ENDIF
      ENDIF
!----------------------------------------------------------------------
!  ARGUMENT IS POSITIVE
!----------------------------------------------------------------------
      IF(Y.LT.EPS)THEN
!----------------------------------------------------------------------
!  ARGUMENT .LT. EPS
!----------------------------------------------------------------------
        IF(Y.GE.XMININ)THEN
          RES=ONE/Y
        ELSE
          RES=XINF
          GOTO 900
        ENDIF
      ELSEIF(Y.LT.TWELVE)THEN
        Y1=Y
        IF(Y.LT.ONE)THEN
!----------------------------------------------------------------------
!  0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
          Z=Y
          Y=Y+ONE
        ELSE
!----------------------------------------------------------------------
!  1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
!----------------------------------------------------------------------
          N=INT(Y)-1
          Y=Y-CONV(N)
          Z=Y-ONE
        ENDIF
!----------------------------------------------------------------------
!  EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
!----------------------------------------------------------------------
        XNUM=ZERO
        XDEN=ONE
        DO 260 I=1,8
          XNUM=(XNUM+P(I))*Z
          XDEN=XDEN*Z+Q(I)
  260   CONTINUE
        RES=XNUM/XDEN+ONE
        IF(Y1.LT.Y)THEN
!----------------------------------------------------------------------
!  ADJUST RESULT FOR CASE  0.0 .LT. ARGUMENT .LT. 1.0
!----------------------------------------------------------------------
          RES=RES/Y1
        ELSEIF(Y1.GT.Y)THEN
!----------------------------------------------------------------------
!  ADJUST RESULT FOR CASE  2.0 .LT. ARGUMENT .LT. 12.0
!----------------------------------------------------------------------
          DO 290 I=1,N
            RES=RES*Y
            Y=Y+ONE
  290     CONTINUE
        ENDIF
      ELSE
!----------------------------------------------------------------------
!  EVALUATE FOR ARGUMENT .GE. 12.0,
!----------------------------------------------------------------------
        IF(Y.LE.XBIG)THEN
          YSQ=Y*Y
          SUM=C(7)
          DO 350 I=1,6
            SUM=SUM/YSQ+C(I)
  350     CONTINUE
          SUM=SUM/Y-Y+SQRTPI
          SUM=SUM+(Y-HALF)*LOG(Y)
          RES=EXP(SUM)
        ELSE
          RES=XINF
          GOTO 900
        ENDIF
      ENDIF
!----------------------------------------------------------------------
!  FINAL ADJUSTMENTS AND RETURN
!----------------------------------------------------------------------
      IF(PARITY)RES=-RES
      IF(FACT.NE.ONE)RES=FACT/RES
  900 GAMMA=RES
!D900 DGAMMA = RES
      RETURN
! ---------- LAST LINE OF GAMMA ----------
      END function gamma


end module cldwat2m_micro
